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SUMMARY 


The general design method for three-dimensional, potential, incompressible or 
subsonic- compressible flow developed in part I of this report is applied to the design of 
simple, unbranched ducts. A computer program, DIN3D1, is developed and five numerical 
examples are presented: a nozzle, two elbows, an S-duct, and the preliminary design of a 
side inlet for turbomachines. The two major inputs to the program are the upstream 
boundary shape and the lateral velocity distribution on the duct wall. As a result of these 
inputs, boundary conditions are overprescribed and the problem is ill posed. However, it 
appears that there are degrees of "compatibility" between these two major inputs and 
that, for reasonably compatible inputs, satisfactory solutions can be obtained. By not 
prescribing the shape of the upstream boundary, the problem presumably becomes well 
posed, but it is not clear how to formulate a practical design method under this circum- 
stance. Nor does it appear desirable, because the designer usually needs to retain control 
over the upstream (or downstream) boundary shape. 

The problem is further complicated by the fact that, unlike the two-dimensional 
case, and irrespective of the upstream boundary shape, some prescribed lateral velocity 
distributions do not have proper solutions. 


1.0 INTRODUCTION 

In part I of this report (ref. 1), a general design method is developed for three- 
dimensional, potential, incompressible or subsonic-compressible flow fields with arbitrary, 
prescribed velocity distributions as a function of arc length along streamlines on the 
boundary of the field. For (the present) part II of this report, a computer program, 
DIN3D1, has been developed for the design of simple, unbranched ducts with uniform 
velocities at the upstream and downstream boundaries and with arbitrary, prescribed 
velocity distributions along streamlines on the lateral boundaries. 

The design of flow fields with satisfactory velocities along the boundary is impor- 
tant for the following reasons: 

(1) Boundary-layer separation losses can be avoided by prescribing velocity distri- 
butions in the direction of flow, along the surfaces of the boundary, that do not decrease 
too rapidly. 

(2) Shock losses in compressible flow and cavitation in incompressible flow can be 
avoided by prescribing velocities that do not exceed certain maximum values dictated by 
these phenomena. 

(3) For compressible flow in ducts, the desired flow rate can be assured by pre- 
scribing velocities that do not result in premature choked flow. 
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However, the first objective of fluid dynamic design is to determine the shape of the 
flow-field boundary for which losses are minimum. For both incompressible and shock- 
free compressible flow, the fluid losses originate at the material surfaces along the flow- 
field boundary, and the magnitude of these losses depends on the velocity distribution 
along these surfaces. The characteristics of a desirable velocity distribution are rela- 
tively well known from boundary-layer theory. 

The main program, DIN3D1, together with 21 major subroutines is described herein. 
The program input and output are also described and several numerical examples are 
presented. 


2.0 SYMBOLS 


All quantities are nondimensional; velocity is expressed as a ratio of the upstream 
velocity; linear quantities, expressed in any consistent unit for input, are made dimen- 
sionless in the program by dividing by the square root of the upstream boundary area; 
names for variables and parameters used in the computer program are not listed. 


A 


B 

C C ,C 0 , 

C i,...,C^ 

c*,c£, 


local continuity parameter, eq. (3.3.2) 

distance between adjacent nodal points of finite-difference star, figs, in 
section 3.2 

local continuity parameter, eq. (3.3.3) 

coefficients in governing finite-difference eqs. (3.3.4) and (3.3.5) 
coefficients given by eqs. (3.4.3) to (3.4.5) 


e 

®1 

®2 

®3 


I, J»K 

i,j,k 

JX,KX 

k 

m 


unit vector 

unit vector in direction of q along streamlines, which are intersections of 
and ti stream surfaces, fig. in section 3.1 

unit vector tangent to intersection of ri stream surface and (p potential 
surface, fig. in section 3.1 

unit vector tangent to intersection of rjf stream surface and <p potential 
surface, fig. in section 3.1 

indices in <p, rjr, and ri directions, respectively 

unit vectors in x, y, and z directions, respectively, fig. in section 3.1 

indices for location of primary streamline, fig. in section 3.5 

adjustment factors for constraints, eqs. (4.2.2) and (4.2.9) 

path length in direction of ¥3 along intersection of stream surface and 
cp potential surface, fig. in section 3.1 



n 


path length in direction of &2 along intersection of tj stream surface and 
<p potential surface, fig. in section 3.1 

Py path length along perimeter of upstream boundary, fig. in section 3.5 

Q In q 

q velocity, expressed as ratio of upstream velocity 

q velocity vector, e^q, fig. in section 3.1 

91 residual error, eq. (3.3.4) 

r radius from axis of duct turn, fig. in appendix A 

s path length along streamline in direction of e^, fig. in section 3.1 

x,y,z Cartesian coordinates in physical space 

a,fi,Y angles of direction cosines in x,y,z space, fig. in section 3.1 

A finite increment 

0 angle with which tJt and T| stream surfaces intersect on potential surface, 

fig. in section 3.1; any angle 

p local static density of fluid as ratio of upstream static density 

<p velocity potential, eq. (3.0.1) 

cp,x|;,Ti curvilinear coordinates in physical x,y,z space, or orthogonal coordinates of 

transformed <p,ty,T] space 

(•) dot-product operator of two vectors 

Subscripts: 

amp amplitude 

D downstream boundary 

max maximum 

min minimum 

o outer radius 

P principal streamline, IP on first fig. in section 7.1 

U upstream boundary 

0,1,.. .,6 grid points in finite-difference star, figs, in section 3.2 

variables or components of variables associated with directions of ei, e~2> 
and ejj, respectively 


1,2,3 
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3.0 GOVERNING DIFFERENTIAL EQUATION AND CONSTRUCTION OF FLOW FIELD 


A major difficulty faced by all classical design methods is that they are boundary- 
value problems in which the velocity is specified along physical boundaries, the shapes of 
which are not known until the problem is solved. In this report, the difficulty was avoided 
by solving the problem in transformed cp,ty,T] space, where <p and ijr are the velocity 
potential and a stream function, respectively, and ti is a second stream function asso- 
ciated with continuity in three-dimensional flow. The velocity distribution can then be 
expressed as a function of <p along the boundary streamlines (lines of constant ip- and ti) 
from the relation 


cp = y*q(s) ds + const (3.0.1) 

where q(s) is the prescribed velocity as a function of arc length s along the boundary 
streamline in physical x,y,z space. 

3.1 Physical x,v,z space . - The flow field at a point in physical x,y,z space has 
two stream surfaces of constant ^ and tj, respectively, that intersect the potential 
surface at 90° and intersect one another at an angle 9 measured on the potential sur- 
face (ref. 1). The directions of these three intersections are given by the unit vectors 
ei, e^, and e^, each defined by its direction cosines cos a, cos E, and cos y. Differ- 
ential lengths along the intersections are given by ds, dm, and dn, as shown in the fol- 
lowing figure. 



cp potential surface-' 




The velocity vector q (e^q) is tangent to the intersection of the and ri stream 
surfaces so that q is normal to the potential surface cp and rp and ti are constant 
along the streamline. 

3.2 Transformed (p,rp,n space . - For a duct, the flow field in transformed <p,tp,T] 
space becomes a cylinder with a cross section 



the same shape as the upstream boundary configuration in x, y,z space, provided that the 
stream surfaces tp and T) at the upstream boundary are defined by lines of constant y 
and z, respectively (ref. 1). Lines of constant ip and ti (paired values) on the lateral 
boundary are streamlines, and the velocity vector q is everywhere parallel to the (p 
axis. The rectangular grid resulting from the intersections of surfaces of constant <p, Tp, 
and T), for various specified values of grid spacing a^, ag, and ag, respectively, is used to 
solve the governing differential equation by finite-difference methods. (The procedure is 
outlined in ref. 1.) Thus, for every internal grid point (numbered 0) in the flow field, at 
which points the finite-difference form of the governing equation must be satisfied, a 
finite-difference star is formed with six adjacent grid points numbered 1 to 6 and spaced 
aj,..., a^ distance away. 
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3.3 Governing differential equation . - From page 28 of part I, the governing, 
second-order, partial-differential equation for the distribution of In q in transformed 
cp,ijr,T) space is 


9 2 ln q 
9cp 2 


9 2 In p 

a 2 
9<p 


9 In sin 9 
9cp 2 


K K, 

n ilr 

2 ' 2 

<1 q 


+ B‘ 


4 A 


/ 9 In B 

aja.a 

+ 9. 2 1 n q\ 

9 In B / 

^ Brlr 

9^ 

at 2 j 

" a<p \ 

/ 9 In A 

9 In q 

+ 9 2 In q\ 

9 In A 

l an 

9ri 

3t? J 

9(p 


SJg-fl 

9<p 


9 In q 
9(p 


9 In B 
9cp 


9 In A 


9(p 


= 0 


(3.3.1) 


where q is the local velocity expressed as a ratio of the upstream velocity (q = 1.0 at the 
upstream boundary), p is the density expressed as a ratio of the upstream static density 
(p = 1.0 at the upstream boundary), 9 is the distortion angle, which is the angle of 
intersection between the ty- and ti stream surfaces, as shown by the figure in section 3.1 
(at the upstream boundary where the specified grid is rectangular, 9 equals 90°, i.e., 


there is no "distortion"), 
and ti stream surfaces (K^ 


and 

and 


K. 

K t 


ti 


are the total curvatures in x,y,z space of the tJt 
equal 0 at the upstream boundary), and A and B 


are the "continuity" parameters defined by (ref. 1) 

a dn • ^ 

A = P dijr sm 9 


(3.3.2) 


dm . 

IS.p^sme 


(3.3.3) 

(A and B equal 1.0 at the upstream boundary, because there d^ = dn and dr\ = dm.) 


In finite-difference form, equation (3.3.1) becomes (ref. 1) 


C c * C 1 Q 1 * C 2 Q 2 + C 3 Q 3 * C 4 Q 4 + C S Q S * C 6 Q 6 ’ C 0 Q 0 


(3.3.4) 


where 






/ 9 In B ^5 

\ ^ a 2 a 5 



+ A 


2/ 8_ln_A 
^ 3ti 




/ 9 In B + d In A\ ( a 4 a l\ 

V 0<p + 0( P )\ a i a 4 ) 
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where Q is In q, the numerical subscripts refer to the six adjacent grid points in the 
finite-difference star shown in the second figure of section 3.2, and the residual error 
0t equals 0 when the governing equation (3.3.1) is satisfied locally. With the coefficients 
Cq, Ci,...,C 6 , Cq at- each internal grid point held constant, equation (3.3.4) is solved 
globally by changing the values of Qq according to standard relaxation procedures. 

These procedures involve repeated passes (IT counter) through the entire flow field, 
starting at the upstream boundary, until the maximum value of 0t in the entire field is 
less than the input value of EPSR. This set of calculations, involving fixed values of the 
coefficients, constitutes one major iteration (ITER counter). The coefficients are then 
recomputed from the new values of Qq, and the procedure is repeated until the solution is 
complete. 


3.4 Direction cosines . - With the velocity gradients known from the solution of the 
governing equation (3.3.1), the nine direction cosines associated with the three unit vec- 
tors e^, and F 3 (fig. in section 3.1) are obtained from their gradients in the (p 
direction starting from the upstream boundary. (The boundary is assumed to lie on a y,z 
plane so that cos ctq, cos 63 , an< * cos Y 3 are 1*0 and. the other six direction cosines are 
zero.) For example, the gradients of the direction cosines of the unit vector e^, which is 
in the direction of the velocity vector "q, are obtained (ref. 1 ) from the components of the 
irrotationality equation normal to the and r| surfaces and from the direction-cosine 
law 


2 2 
cos a + cos 


B + cos Y = 1 


(3.4.1) 


Thus, 
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and 



(3.4.3) 



(3.4.4) 


Equations for the gradients of the direction cosines of &2 an< * e 3 the <P direc- 
tion are obtained in part I of this report (ref. 1) in a similar manner. In (the present) part II, 
however, the expressions for C3 and C4 (eqs. (22f) and (23f) in part I) have been 
reformulated, by making use of the continuity equations (16c) and (16d) in reference 1, to 
give 



Finally, again by making use of the continuity and irrotationality equations in part I, 
equations can be developed (section 5.2) for the gradients of the direction cosines of the 
unit vector e^ in the and T| directions on potential surfaces. These gradients are 
used in subroutine ANGL (section 5.5). 


3.5 Construction of flow field in x,v,z space . - With the velocity distribution 
known from the solution of equation (3.3.1) in <p,*|r,Ti space, and with the distribution of 
the direction cosines likewise known from section 3.4, the shape of the flow field in x,y,z 
space can be constructed. The boundary of this flow field constitutes the design of the 
duct. 


The construction starts in x,y,z space with the arbitrarily specified shape of the 
upstream boundary on the yu>zu plane at xu = 0, where, for nondimensional variables (as 
defined in ref. 1), r|r and ri equal yu and zjj> respectively, because the grid is rectan- 
gular (but not necessarily square). 
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From each intersection of the grid lines in the figure, a streamline (with constant 
paired values of and t}) extends to the downstream boundary. The x,y,z coordinates 
of the streamline at each successive potential surface (constant cp) are obtained by 
integrating the following equations (eqs. (26a), (26b), and (26c), ref. 1): 


and 



(3.5.1) 


(3.5.2) 


(3.5.3) 


3.6 Alternative construction of flow field in x,v,z space . - An alternative method 
for constructing the flow field is to select one streamline (designated by the indices JX 
and KX as shown on the figure in section 3.5), obtained from equations (3.5.1) to (3.5.3), 
and to use this primary streamline as a backbone from which to obtain the x,y,z coor- 
dinates of every grid point on each successive potential surface by integrating the fol- 
lowing equations (ref. 1) in the ip direction on a potential surface, 



and in the ti direction on a potential surface, 



(3.6.1) 


(3.6.2) 


(3.6.3) 


(3.6.4) 


(3.6.5) 
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and 



The continuity parameters A and B in equations (3.6.1) to (3.6.6) are computed by the 
program from equations (3.3.2) and (3.3.3), respectively, by using values of An and Am 
from the previous iteration. 

The location JX,KX of the primary streamline is arbitrary. However, results 
should be best for locations near the center of gravity of the upstream boundary and 
should be biased somewhat toward the boundary streamlines with the highest prescribed 
velocities if the duct bends. For solutions with planar symmetry, the computer program 
DEN3D1 requires that JX,KX be on the plane of symmetry. 

In the program, both methods (sections 3.5 and 3.6) are used to find the x, y,z 
coordinates. This is further discussed in section 4.7, where the input coefficient CAVP 
is introduced to allow a weighted average of the two methods. The first method 
(section 3.5) is used in subroutine VARI (section 5.2), where CAVP also is used, and the 
second method (section 3.6) appears in subroutine POTS (section 5.16). 


4.0 ILL-POSED NATURE OF DESIGN METHOD WHEN APPLIED TO DUCTS 

The design method applied to ducts requires two major inputs: (1) the upstream 
boundary configuration and (2) the velocity distribution on the lateral boundary. The 
lengths As of all streamlines on the boundary are precisely fixed because along each 
streamline 



(4.0.1) 


where q is a known function of (p from equation (3.0.1) or, alternatively, is specified 
directly as a function of (p. Thus, for various upstream boundary configurations, which 
for a uniform (constant) upstream velocity with parallel flow must be plane, it appears 
unlikely that the downstream potential surface can also be plane with parallel streamlines 
normal to the surface, as required by the design method. If this is the case, boundary 
conditions are overprescribed and the design problem is ill posed. 

For every prescribed upstream boundary configuration that lies on a flat, potential 
surface as assumed by the design method (section 7.2), there is an infinity of compatible 
velocity distributions that could exist on the lateral boundary, because there is an infinity 
of lateral boundary configurations for any upstream boundary shape. However, this con- 
sideration does not rule out the possibility of an infinity of lateral velocity distributions 
that are not compatible with the prescribed upstream boundary configuration. 

Now, for a given duct configuration (completely specified) with upstream and down- 
stream regions extended so that the upstream and downstream boundaries are flat poten- 
tial surfaces, as assumed by the design method, a specific velocity distribution exists 
throughout the flow field and in particular on the lateral boundary. Presumably, this 
velocity distribution is unique to this duct shape (e.g., pp. 14 to 41, ref. 2); it then follows 
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that for a given lateral velocity distribution there is a unique upstream boundary config- 
uration. Thus, although for a given upstream boundary an infinity of lateral velocity 
distributions exists, as discussed in the previous paragraph, for a given lateral velocity 
distribution there is only one compatible upstream area configuration. It is concluded 
that the general design method when applied to ducts is ill posed because the boundary 
conditions are overprescribed. However, fortunately, there are "degrees of incompati- 
bility," as considered in the next section. 

4.1 Compatibility between prescribed upstream boundary configuration and lateral 
velocity distribution . - For a given lateral velocity distribution, some upstream boundary 
configurations are less compatible than others. For example, a stellated upstream con- 
figuration such as 



would be "highly incompatible" with a lateral velocity distribution corresponding to the 
flow through an elbow of constant, circular cross section. An elliptical upstream con- 
figuration with moderate aspect ratio should be "highly compatible" with such a lateral 
velocity distribution. 

Program DIN3D1 has been so constructed that stream-tube areas are adjusted to 
local velocities (by means of the continuity parameters A and B; section 4.2 (con- 
straint 6)), so that, unless the upstream boundary configuration is "highly incompatible," 
the continuity condition is essentially satisfied (the downstream- are a error is 
printed) for each of a relatively large number of major iterations (ITER; end of 
section 3.3). The solution at first converges for each successive major iteration (as evi- 
denced by decreasing maximum Si in the flow field). Eventually, because boundary 
conditions are overspecified, the solution must diverge and fail. For "highly incompatible" 
upstream boundary configurations, divergence is rapid and occurs after only a few major 
(ITER) iterations. For "highly compatible" cases, divergence is gradual and occurs only 
after many iterations. Thus, excellent approximate solutions are obtained by stopping the 
calculations before, or shortly after, divergence begins (section 4.4). 

Finally, the program includes options (IVEL equals 2 or 3) for lateral velocity dis- 
tributions that tend toward compatibility with the prescribed upstream boundary config- 
uration. These "equilibrium" velocity distributions are based on the velocity distribution 
in ducts of constant cross section and very large turning angle. Under these conditions, 
near the middle of the turn, the velocity distribution on potential surfaces becomes a free 
vortex (qr = constant) with r measured from the axis of the turn. A lateral velocity 
distribution around the periphery of each potential surface, based on this free-vortex 
distribution, and with r related to the upstream surface configuration, constitutes the 
"equilibrium" velocity distribution on the lateral boundary (appendix A). 
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4.2 Constraints on calculation procedure . - Because of the ill-posed nature of the 
method when applied to ducts, it is beneficial to guide the calculations by imposing the 
following six constraints, all of which would be satisfied automatically in a well-posed 
case: 


(1) The x,y,z coordinates of every internal grid point in the flow field are computed 
by two methods (sections 3.5 and 3.6), and the results are averaged according to the input 
value of CAVP, the decimal fraction of the first method that enters into the weighted 
average. Thus, 


0.0 < CAVP <1.0 (4.2.1) 

Although the optimum value of CAVP probably varies with the complexity of the 
upstream boundary configuration and with the lateral velocity distribution, a value of 0.5 
is generally satisfactory (also see section 4.6). 

(2) During the iterative calculations, values of the direction cosines computed from 
their gradients (eqs. (3.4.2), e.g.) can become greater than 1.0. If this occurs, the value is 
set equal to 1.0 by the program. 

(3) Also, during the iterative calculations, the sum of the squares of the direction 
cosines may not equal 1.0, as required by equation (3.4.1). When this occurs, the program 
changes each cosine value by a factor k, where 

k 1 (4.2.2) 

/ o 2 2 

I cos a + cos R + cos yl 

(4) At every interior grid point, the unit vector e^, which is tangent to the 
intersection of the T| stream surface and the velocity potential surface (p, must be 
normal to the unit vector e^, which is in the direction of the velocity (fig. in 
section 3.1). Thus, 


*1 ' *2 -° 


(4.2.3) 


or 


cos a cos a + cos J3 cos J3 + cos y cosy =0 (4.2.4) 

1 L l 2* i z 

When this relation is not so, the direction cosines of &2 are changed by the program so 
that &2 becomes normal to e^ and the plane of e^ and &2 remains unchanged. This 
same "normality" condition is imposed on the direction cosines of the unit vector F3 in 
the same figure. (Also, see appendix B.) 

(5) Also, from the figure in section 3.1, 

^2 ‘ cos 9 (4.2.5) 


from which 


cos 9 = cos cos + cos R^ cos B^+ cos y^ cos y^ 


(4.2.6) 
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From trie following figure, which shows a stream tube bounded by adjacent surfaces of 



constant t|t and tj, the value of 9 must be greater than 0.0° and less than 180.0°; 
otherwise the stream- tube area becomes zero or negative. Thus, if from equation (4.2.6) 
the absolute value of cos 0 is greater than 0.9962, the program sets cos © equal to 
+0.9962 so that the "distortion" angle 9 lies in the range 

5° <©<175° (4.2.7) 

(6) Finally, the values of the continuity parameters A and B, as computed by 
equations (3.3.2) and (3.3.3), respectively, are changed by the same factor k so that the 
following continuity condition (eq. (lOd), ref. 1) is satisfied: 



(4.2.8) 


from which 



(4.2.9) 


4.3 Mode of failure . - For those cases where total failure of the solution is 
approached after a sufficiently large number of major iterations (ITER), which number 
depends on the compatibility condition discussed in a previous section, 4.1, this failure 
usually occurs in the downstream region. Here the flow field in x,y,z space distorts as 
shown by the following examples. 
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The lateral surfaces diverge, and the "distortion" angle 9 (constraint (5), section 4.2) may 
vary rapidly and greatly from its initially undistorted value of 90° at the upstream bound- 
ary (fig. in section 3.5). All of these distortions result from an accumulation of unreal 
values of In q near the downstream boundary. This accumulation appears to result from 
the relaxation procedure, which always starts at the upstream boundary and marches 
through to the exit, continually pushing the effects of the ill-posed problem toward the 
downstream boundary. The constraints discussed in section 4.2 maintain an apparently 
well-behaved flow field elsewhere. (In section 4.8 the distortion near the downstream 
boundary is found to be essentially independent of the extent of the downstream region, 
supporting the preceding reasoning.) 

4.4 Input option ISOLV . - As stated in section 4.1, excellent approximate solutions 
can be obtained for prescribed upstream boundary configurations and lateral velocity 
distributions that are moderately compatible. These solutions are achieved by stopping 
the calculations before or shortly after divergence begins. Input option ISOLV = 1 
assumes that, because of the constraints discussed in section 4.2, an adequately converged 
solution is achieved after four or more major iterations (ITER), provided that the error in 
exit flow area, expressed as a decimal fraction, is less than 0.0033. Also, for ISOLV = 1, 
if these criteria are not met, the calculations are then stopped and the solution is printed 
out when the computed value of exit-area error changes sign - provided that the value of 
ITER is greater than 8. In this latter case, the solution may not be acceptable if the 
exit-area error (intermediate printout) is changing rapidly as the result of impending 


16 



failure. If none of the criteria are met, the solution continues until the number of major 
iterations equals the input value of ITERMX or until the solution fails entirely. 

For an input value of ISOLV = 0, the solution continues until (1) the number of major 
iterations (ITER) equals the input value of ITERMX, (2) provided that ITER > 4, the max- 
imum value of 0t (eq. (3.3.4)) in the entire flow field at the beginning of a major iteration 
is less than or equal to 0.0020, or (3) the solution fails. 

4.5 Effect of ITER on solution . - Because of the ill-posed nature of the design 
method when applied to ducts, all examples must eventually fail (with rare exceptions in 
simple cases) as the number of major iterations (ITER) increases indefinitely. Thus, as 
shown, the solution changes with ITER. The changes are most pronounced in the 



ITER = 4 ITER - 12 ITER - 24 



ITER - 48 



ITER * 72 


downstream- area configuration and to a lesser degree in the turning angle of the duct. 

It is suggested that, because the exit velocity is normal to the downstream boundary and 
therefore not influenced by its shape, large changes in the exit-area configuration can 
result from only minor changes in the lateral velocity distribution (appendix C). Likewise, 
for the same lateral velocity distribution during an approximate solution, small changes in 
the. downstream boundary shape (but not in its area) can be expected from one major 
iteration (ITER) to the next. Because for the larger values of ITER the solution is 
approaching failure, the lower values of ITER are believed to give better approximate 
solutions. For ISOLV = 1, the solution is stopped at ITER = 4, provided certain criteria 
are met (section 4.4). 
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4.6 Effect of CAVP on solution . - In sections 3.5 and 3.6, two methods are 
discussed for computing the coordinates of the flow field in physical x,y,z space. The 
first method (section 3.5), which results in correct streamline lengths, determines x, y, 
and z by integrating along streamlines between adjacent potential surfaces. (For this 
method, the streamline lengths are correct, but the continuity condition may not be 
satisfied because of the ill-posed nature of the problem and its forced solution.) The 
second method (section 3.6), for which continuity is satisfied, starts from the "primary" 
streamline location (JX,KX; fig. in section 3.5) and determines the coordinates by 
integrating in the t|r and ri directions on potential surfaces. (For this method, the 
continuity condition is satisfied by the essentially correct values of the continuity 
parameters A and B, but the streamline lengths may not be correct.) When moving from 
one potential surface to the next, these two sets of x,y,z coordinates are averaged by 
the input value of CAVP, which is the decimal fraction of the first set that enters into 
the weighted average of the two sets. As shown in the following figures, the effect of 
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CAVP on the solution for ITER = 72 is to go from the case where streamlines are normal 
to the downstream potential surface but not parallel with each other (CAVP = 0.0) to the 
case where the streamlines are parallel but not normal to the potential surface (CAVP = 
1.0). A clearer picture of the difference is given by two side views: 




For less extreme examples (ITER << 72), the two methods of computing the x,y,z 
coordinates should give more nearly equal results, and a value of 0.5 for CAVP should 
usually be satisfactory. For more complex upstream boundary configurations and for 
problems involving less compatibility between the upstream boundary shape and the pre- 


19 


scribed lateral velocity distribution, where the solution starts to diverge at relatively low 
values of ITER, a better input value for CAVP might be as low as 0.2. However, in 
some cases, because of the numerical integration procedure used by the second method 
(section 3.6) on the potential surface, the duct wall may develop slightly rippled regions. 
The ripples can be reduced by increasing the input value of CAVP or eliminated by 
setting CAVP equal to 1.0. 

4.7 Effect of duct turning angle on solution . - As might be expected, the greater 
the duct turning angle A0, the lower the value of ITER at which the solution begins to 
diverge. Or, as a corollary, for the same value of ITER, other things being equal, the 
greater A0, the greater the distortion (if any) near the downstream boundary, as shown 
by the following figures: 



A0 ~ 420° 



A0 ~ 720° 
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4.8 Effect of extent of downstream region on solution . - The downstream region is 
that part of the duct near the downstream boundary where the lateral velocity is constant 
and equal to its value at the downstream boundary. Because solutions apparently start to 
diverge and distort near the downstream boundary, a question arises as to whether the 
extent of the downstream region influences the magnitude and type of this distortion. As 
the following figures indicate, the distortion is apparently not influenced appreciably by 
the extent of the downstream region. This observation supports the discussion in 
section 4.3 regarding the mode of failure, in program DIN3D1, for ill-posed problems. 



4.9 Effect of complexity of upstream boundary configuration on solution . - Exam- 
ples of solutions with simple upstream boundary configurations and others with relatively 
complex configurations are given by the figures on the next page. Solutions for complex 
configurations require special care in prescribing lateral velocity distributions. Simple 
upstream configurations involve no particular difficulty. 

4.10 Remarks . - In concluding this major section on the ill-posed nature of the 
general design method when applied to ducts, it is noted that the method is not ill posed 
when applied to the external design problem (i.e., to the design of bodies with prescribed 
velocities in infinite space). In this case, the velocity on the outer boundaries is every- 
where constant, provided only that the upstream and downstream boundaries, which can 
have any suitable shape (e.g., circular), are sufficiently large. This situation eliminates 
the compatibility problem between prescribed boundary configuration and prescribed 
velocity distribution, but other problems, such as body closure, remain. 

Thus far, it has been implied that the general design method when applied to ducts 
becomes well posed if the upstream boundary configuration is not specified. Under this 
circumstance, it is not clear how, and may not even be possible, to carry out a practical 
design procedure. Nor does it appear desirable, because the designer usually needs to 
retain control over the upstream boundary configuration (which, by reversing the direction 
of flow, becomes the downstream configuration, assuming that configuration needs to be 
controlled). 


21 




The difficulty of the three-dimensional duct design problem is further increased by 
the likelihood that, irrespective of the upstream boundary configuration, and unlike the 
two-dimensional case, proper solutions do not exist for every lateral velocity distribution. 
The situation is considered in more detail in appendix C. 


5.0 BRIEF DESCRIPTION OF COMPUTER PROGRAM DIN3D1 

Program DIN3D1 is written in standard Fortran IV. Double precision is required for 
computers with 32-bit words. In addition to the main program, there are 21 major sub- 
routines, 3 minor subroutines, and 8 external functions. The main program and the major 
subroutines are briefly described in this section. 

5.1 Main program . - The main program governs the solution, as shown on the sim- 
plified flowchart. The integer I designates a potential surface, starting with 1 at the 
upstream boundary and ending with NT at the downstream boundary; the integer IT 
counts the number of passes through the entire flow field (1 = 2 to I = NI-1= Nil), all 
with the same global set of coefficients in the finite-difference equation (3.3.4); and the 
integer ITER counts the number of major iterations, each with an improved set of coef- 
ficients determined from the previous major iteration. 



Simplified flowchart of main program 
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In general, the design procedure in program DIN3D1 is as outlined under Numerical 
Procedure on pages 56 and 57 in part I (ref. 1) of this report. The overall approach is to 
solve (in subroutine RELAX) the finite-difference form of the governing differential 
equation (3.3.4) everywhere on one potential surface at a time, starting at I = 2 (the first 
surface downstream from the upstream boundary) and marching through the entire flow 
field to I = Nil. This procedure is continued (counter IT), with the same global set of 
coefficients in the finite-difference equation (3.3.4), until the maximum residual error 
RESMAX (flowchart) encountered anywhere in the flow field is less than the input value 
of EPSR, or until IT equals the input value of ITMAX. At this point, a new set of 
coefficients is generated (in subroutine COEF, by using major parameters determined in 
subroutine VARI); the error (ERRAR) in the downstream boundary area (expressed as a 
decimal fraction of the correct value) is computed (flowchart); and the procedure is 
repeated. This process is continued (counter ITER) until the maximum residual error 
REMAX (flowchart) encountered anywhere in the flow field on the first pass with a new 
set of coefficients is less than 0.002 (provided ITER > 4), or until the value of ITER is 
equal to the input value of ITERMX (flowchart). This concludes the solution of the 
governing equation (3.3.4), which, as indicated on the flowchart, may also be concluded 
sooner if the input value of ISOLV is equal to 1 (section 4.4). 

The flowchart for the main program involves 10 (of 21) major external subroutines. 

These are described shortly. 

5.1.1 Input ICONX. - In the main program, to achieve better accuracy and to speed 
up the solution, the new global set of coefficients (for the governing equation) calculated 
after each major iteration (ITER) is iterated (counter ICON). This iteration is continued 
until ICON is equal to the input value of ICONX (flowchart). However, to further 
shorten the solution time, the value of ICONX is reduced by 1 after every seven major 
iterations (ITER) until a minimum value of 3 is attained, after which ICONX remains 
constant. (If the input value of ICONX is less than 3, it is changed to 3 in the main 
program and remains constant.) 

5.2 Subroutine VARI . - Except for the main program, subroutine VARI is the most 
important routine in program DIN3D1. After each major iteration (ITER counter in main 
program), this subroutine determines, at every grid point in the flow field, the direction 
cosines of the unit vectors e"i, e"2, and F3 (section 3.4); the x, y, and z coordinates (section 
3.5); and the parameters A, B, and 9 and the coefficient Cc (sections 3.3 and 4.2 
(constraint 5)). The procedure is outlined in the simplified flowchart. In general, this 
procedure is as outlined in part I of this report (pp. 56 and 57, ref. 1). For a given poten- 
tial surface I, the overall approach in subroutine VARI is to compute the direction cosines 
of the unit vectors e^, l^, and eg on potential surface I + 1 from their known values on 
potential surface I and from their derivatives with respect to cp on both the I and I + 1 
surfaces, as given by equation (3.4.2), for example. Because the derivatives of the direc- 
tion cosines on surface I + 1 depend on the direction cosines themselves, the procedure is 
iterated three times (IT2.EQ.3, flowchart). Afterward (IT1.GT.1, flowchart) the direction 
cosines of the unit vector e^ are determined in subroutine ANGL (flowchart) by a new 
method based on their known values at the primary streamline location (specified by the 
input values of JX and KX; section 3.6 and fig. in section 3.5) on potential surface I + 1 
and from their derivatives with respect to and tj. 
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Simplified flowchart of subroutine VARI 


These ^ and ti derivatives of the three direction cosines of ¥]_ are obtained by the 
following procedure. The derivatives with respect to r|r are obtained from the simulta- 
neous solution of the following three equations: (1) the T|r derivative of the direction- 
cosine law (eq. (3.4.1)), (2) the continuity equation (16c) from reference 1, and (3) the 
equation (5.2.1). This equation is obtained by adding the irrotationality equations (14c) 
and (14d) from reference 1 to obtain 3 cos 0/3cp, after which equations (13g), (16c), and 
(16d), also from reference 1, are introduced to give 


e 


3 


3e. 


1 1 
3ty - 2B 


3 cos Q 
3cp 


- cos 0 


9 In q 
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(5.2.1) 


25 
















In a similar fashion, the derivatives with respect to ti are obtained from (1) the 
direction- cosine law, (2) the continuity equation (16d), and (3) equation (5.2.1) combined 
with equation (13g) from reference 1 to give 


e 2 


d e j ^ _l_ 9 cos 6 

9n = 2A 9(p 


- cos 9 


( 0 3 In q 

9 In A 

9 In B\' 

(, 2 3<P 

+ 9<p + 

9q> )_ 


(5.2.2) 


This new method for finding the direction cosines of e^ is needed to achieve truly 
parallel flow in the downstream boundary region. Using this method in subroutine ANGL, 
the procedure (in subroutine VARI) is to iterate NTRY additional times (where, if the 
input value of NTRY is less than 2, the program sets NTRY equal to 2). 


Finally, as shown in the flowchart, subroutine VARI determines the x,y,z coor- 
dinates of the physical flow field (sections 3.5 and 3.6), the distortion angle 9 (section 4.2 
(constraint 5)), the continuity parameters A and B (section 4.2 (constraint 6)), and the 
coefficient Cq, which is required by the finite-difference equation (3.3.4). 

Subroutine VARI is called from the main program (section 5.1) and from subroutine 
PUTOUT (section 5.18). 

5.2.1 ICX, ICY, and ITH counters. - If the three simultaneous equations (e.g., 
eq. (3.4.2)) for the derivatives of the three direction cosines of each of the unit vectors 
*i, e 2 , ° r 63 are not independent, their determinant D will become zero. Thus, in 
subroutine VARI, if -0.00001 < D < 0.00001, D is set equal to +0.00001, and the counter 
ICX is increased by 1. (This condition can occur when the solution is diverging, and 
failure usually occurs soon after.) If, for a given potential surface I, the value of ICX is 
greater than zero, a CAUTION note appears in the intermediate printout, or if the value 
is greater than 10, another note appears and the solution is stopped. 

If the absolute value of any direction cosine, obtained in subroutine VARI from the 
derivatives of the three direction cosines, is greater than 1.0, that value is changed to 
+ 1.0, and the counter ICY is increased by 1. If, for a given potential surface I, the value 
of ICY is greater than zero, a CAUTION note appears in the intermediate printout. If 
the value is greater than 10, another note appears and the solution is stopped. 

Also, in subroutine VARI, the "distortion" angle 9 is constrained to values between 
5° and 175° (section 4.2 (constraint 5)). If 0 is less than 5° or greater than 175°, cos 9 is 
set equal to 0.9962, and the counter ITH is increased by 1. If, for a given potential sur- 
face I, the value of ITH is greater than zero, a CAUTION note appears in the interme- 
diate printout. If the value is greater than 10, another note appears and the solution is 
stopped. 

5.2.2 Averaging coefficients CAVD, CAVN, CAVX, CAVY, and CAVZ. - In 
subroutine VARI, during the iterations involving direction cosines and their derivatives, 
the new values are averaged with the previous values by the input values of CAVX and 
CAVD, respectively, where CAVX and CAVD are the decimal fractions of the previous 
(old) values entering into the weighted average. 

Likewise, the input values of CAVN, CAVY, and CAVZ are the decimal fractions 
of the previous values of the continuity parameters A and B, the coefficient Cq 
( eq. (3.3.4)), and the cosine of the "distortion" angle 0, respectively, that enter into the 
weighted average with the respective new values. 
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For design problems involving lateral velocity distributions to achieve the desired 
shapes of the downstream boundary, the input value of CAVN can be as as low as 0.0, 
and certainly no higher than 0.1. In effect, the new values of A and B are not averaged 
with the previous values. 

5.3 Subroutine AERIA . - Given the incremental lengths SI, ..., S6, shown in the 
figure, subroutine AERIA determines the incremental area DARE A of an incremental, 



curved surface bounded by four, essentially straight, incremental lines. The area is 
divided into two sets of triangles by S5 and S6. For each of these four triangles 


AA 


= [s (s - aXs - b)(s - c)] 


1/2 


(5.3.1) 


where 


a + b + c 
S = 2 


(5.3.2) 


Thus, s is the semiperimeter of the triangle, and a, b, and c are the lengths of its three 
sides. Subroutine AERIA is called from subroutine ERIA, where it is used to compute the 
area of potential surfaces. Starting with subroutine AERIA, the subroutines are discussed 
in alphabetical order. 

5.4 Subroutine AKA . - Subroutine AKA assures that the sum of the direction 
cosines squared is equal to 1.0. (section 4.2 (constraint 3)). Subroutine AKA is called from 
subroutines ANGL, FINIS, and VARI. 


5.5 Subroutine ANGL . - On potential surface I + 1, subroutine ANGL determines 
the distribution of the direction cosines of "ej, starting from the location of the primary 
streamline (input values of JX and KX) and integrating along lines of constant t|r and ri 
(eqs. (5.2.1) and (5.2.2)). Subroutine ANGL is called from subroutine VARI. 

5.6 Subroutine BOUND . - The physical x,y,z coordinates of the flow field at all 
interior grid points are determined by subroutine VARI. Using these values on a given 
potential surface I, subroutine BOUND extrapolates to determine the coordinates XB, 

YB, and ZB (and the velocity QB) at every contour point along the boundary of the 
potential surface. Subroutine BOUND is called from subroutine PUTOUT. 

5.7 Subroutine COEF . - On potential surface I, subroutine COEF determines the 
values of the coefficients CO, Cl, C2, etc., in the finite -difference form of the governing 
differential equation (3.3.4). Subroutine COEF is called from the main program. 
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5.8 Subroutine ENGL . - On potential surface I + 1, subroutine ENGL determines 
the and ti derivatives of the direction cosines of e^. These derivatives are required 
to obtain the direction cosines of e^. Subroutine ENGL is called from subroutine ANGL. 

5.9 Subroutine ERIA . - On potential surface I, subroutine ERIA computes the flow 
area (dimensional) of the potential surface. Subroutine ERIA is called from subroutine 
FINIS. 

5.10 Subroutine FINIS . - Starting at potential surface NI - 1, subroutine FINIS 
determines the values of A, B, THET, and C£> and of the x,y,z coordinates at every 
internal grid point on the downstream potential surface NI. This subroutine assumes that 
at the downstream boundary the <p derivatives of all direction cosines are zero. Sub- 
routine FINIS obtains the downstream flow area at NI - 1 by calling subroutine ERIA. 
Subroutine FINIS is called from the main program and from subroutine PUTOUT. 

5.11 Subroutine FIRST . - On potential surface I, subroutine FIRST establishes the 
values of the variables appearing in the coefficients (eqs. (3.3.5)) of the finite-difference 
form of the governing differential equation (3.3.5). Subroutine FIRST is called from the 
main program. 

5.12 Subroutine FLAR . - Subroutine FLAR computes the sum (EXFLAR) of all 
incremental flow areas bounded by four internal grid points on potential surface I. Sub- 
routine FLAR is called from subroutine ERIA, which adds to the value of EXFLAR all of 
the incremental areas adjacent to the potential surface contour. 

5.13 Subroutine GRID . - Subroutine GRID determines the area of the upstream 
boundary surface and the distance P around its contour. This subroutine also determines 
the grid spacings a2, ag, and ag (second fig. in section 3.2) on the potential surfaces. 
Further details regarding the potential surface grid are given in section 7.1 of this report. 
Subroutine GRID is called from subroutine PUTIN. 

5.14 Subroutine ONE ST 1 . - On potential surface I, during the first major iteration 
(ITER = 1 only), subroutine ONEST1 determines the parameters A, B, Cq, and THET at 
every internal grid point. This subroutine assumes that the total curvatures and 

of the stream surfaces are zero, that cos 9 (i.e., THET) is also zero, that 

/ n \ 1/2 

A=/J] (5.14.1) 

and that 


B = A 


(5.14.2) 


Subroutine ONEST1 is called from the main program. 

5.15 Subroutine PAR AM . - On potential surface I, subroutine PAR AM determines 
the values of all variables required to compute the coefficient Cq in equation (3.3.4). 
Subroutine PARAM is called from subroutine VARI. 

5.16 Subroutine POTS . - On potential surface I + 1, subroutine POTS determines 
the distribution of the x,y,z coordinates of the internal grid points, starting from the 
known values of x, y, and z for the primary streamline at JX,KX (fig. in section 3.5) and 
integrating in the i|r and ti directions (section 3.6). Subroutine POTS is called from 
subroutine VARI. 
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5.17 Subroutine PUTIN . - Subroutine PUTIN reads and writes all input data required 
by program DIN3D1. Detailed descriptions of these data are given in section 7.3. Sub- 
routine PUTIN is called from the main program. 

5.18 Subroutine PUTOUT . - Subroutine PUTOUT prints those results of the solution 
that are requested by subroutine PUTIN. There are three output tables: (1) output table I 
(internal grid points), (2) output table II (coordinate points along contours of selected 
potential surfaces), and (3) output table III (coordinate points along selected streamlines). 
For details of the data reported in each of these tables, see section 8.0. Subroutine 
PUTOUT is called from the main program. 

5.19 Subroutine RELAX . - At every interior grid point on potential surface I, sub- 
routine RELAX reduces the residual error & (section 3.3) to an absolute value less than 
the current value of EPSX by varying the value of Qg in the finite-difference equation 
(3.3.4). (The initial value of EPSX, which is 400 times the input value of EPSR, 
decreases 12 times by a factor of 0.5, after which its final value is approximately 0.1 
times the input value of EPSR.) This subroutine uses an overrelaxation coefficient, the 
input value of which is O RELAX. The relaxation process is continued until the maximum 
residual error everywhere on the potential surface is less than 0.1 times the input value 
of EPSR or until ITX is equal to the input value of ITXMAX, whichever occurs first. 

During the relaxation process the program marches across the potential surface first 
from left to right (increasing index J), then from top to bottom (decreasing ti index 
K), next from right to left (decreasing J), and then from bottom to top (increasing K). 
Subroutine RELAX is called from the main program. 

5.20 Subroutine RESID . - Subroutine RESID determines the residual error at every 
interior grid point on a given potential surface I. Subroutine RESID is called from the 
main program. 

5.21 Subroutine START . - For ITER greater than 1, subroutine START establishes 
the values of certain variables Cej., ^ > ®3> x > Y> z » Q> RHO, THET, A, and B) at potential 
surface 1 = 1. Subroutine START is called from the main program and from subroutine 
PUTOUT. 

5 .22 Subroutine VELD . - From various input data, subroutine VELD determines the 
input velocity distribution on the lateral boundary of the flow field. From this distribu- 
tion, it estimates initial values of velocity at all interior grid points. For more detailed 
discussion of the input velocity distribution on the lateral boundary, see section 7.2. 
Subroutine VELD is called from subroutine PUTIN. 


6.0 MISCELLANEOUS FEATURES OF PROGRAM 

Various special features of the program, in addition to those already discussed, are 
described in this section. These features relate mainly to user options and to input param- 
eters affecting the running time and accuracy of the calculations. 

6.1 Option IFLUID . - The input option IFLUID relates to the type of fluid used in 
the duct design. At present, provision has been made for two types: incompressible fluids 
(IFLUID = 1), and perfect gases (IFLUID = 2). For incompressible fluids, no additional 
inputs are required (e.g., the fluid density is not required). For compressible fluids, two 
additional inputs are required: the upstream Mach number AMU, and the ratio of specific 
heats GAM. 
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affects the turning angle, and of course, the magnitude of the downstream area. The 
incompressible solution and that for AMU = 0.01 are essentially equal. 

If other types of fluid are required, additions can be made to the code in subroutine 
PUTIN. Of course, appropriate additions to the code must be made throughout the pro- 
gram, wherever the static density ratio p appears. 

6.2 Option ISYM . - Many duct designs have planar symmetry, that is, the duct 
shape on one side of the plane is a mirror image of that on the other. If the prescribed 
lateral velocity distribution has planar symmetry, so also will the resulting design, pro- 
vided that the prescribed upstream boundary configuration is also symmetrical about the 
plane. For cases involving planar symmetry, provision is made in the code for solving only 
one of the two flow fields on either side of the plane of symmetry. This provision cuts the 
running time roughly in half, or alternatively permits a finer grid with the existing 
21 -by- 3 6 array size. 

The input option ISYM relates to three types of symmetry: first (ISYM = 1), there 
is no symmetry, or if planar symmetry exists, it is not made use of; second (ISYM = 2), 
there is, as shown in the following figures, 



ISYM ■ 1 



ISYM * 2 
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symmetry about a plane of constant y (i.e., a plane normal to the y axis of the upstream 
boundary (fig. in section 3.5)); and third (ISYM = 3), there is, as shown in the following 
figure, symmetry about a plane of constant z. 


z 


ISYM =3 



Both types of planar symmetry have been introduced because in program DIN3D1 
the 21-by-36 arrays corresponding to the y and z directions, respectively, are not 
square. 


6.3 Option IP LOT . - Because the output for three-dimensional solutions is usually 
very large, some sort of graphics display is almost a necessity. In program DEN3D1, the 
display data consist of the x,y,z coordinates of selected points along the contour of 
every selected potential surface IPS(200). These same points correspond to selected 
streamlines ISL(200) on the lateral boundary. The result is a three-dimensional plot of 
the lateral boundary of the flow field (i.e., the duct surface) consisting of a network of 
potential surface contours and streamlines as shown by figures in this report. (These same 
data are also printed in output table II (contour data, section 8.3) and output table III 
(streamline data, section 8.4).) 

These display data are obtained by setting the input value of option IPLOT equal 
to 1. (For no graphics display data, the input value of IPLOT is zero.) These graphics 
data, which normally go to tape or disk, are provided for at the end of subroutine PUTOUT 
(section 5.18). It is assumed that a three-dimensional graphics program is available to the 
user, and only the following raw data, in the order presented, are supplied by program 
DIN3D1: 

NI total number of potential surfaces from upstream boundary to 

downstream boundary. (NI is computed by the program from the 
input value of NP (number of potential surfaces for which data 
are specified) and NSD (number of equal subdivisions along the 
principal streamline between each of the NP potential surfaces); 
NI = (NP - 1) NSD + 1.) 

NCP total number of contour points around each potential surface. 

(Each contour point corresponds to a streamline on the lateral 
boundary, so there are NCP streamlines.) 
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NPS 


number of potential surfaces to be plotted and printed out in 
table II (maximum, 200) 


NSL 


number of streamlines to be plotted and printed out in table III 
(maximum, 200 ) 


IPS(200) 

ISL(200) 


index values IY of NPS potential surfaces to be plotted and 
printed (maximum, 200 ) 

index values IX of NSL streamlines to be plotted and printed 
(maximum, 200 ) 


XB(ISL(IX),IPS(IY)) 

YB(ISL(IX),IPS(IY)) 

ZB(ISL(IX),IPS(IY)) 


x,y,z coordinates of contour points around potential surfaces 
and corresponding points along streamlines on lateral boundary 


6.4 Overrelaxation factor (O RELAX) . - At any point in the <p,iJf,T| flow field, the 
residual error St in the governing equation (3.3.4) can be eliminated by an incremental 
change AQq in the local value of Qq. From equation (3.3.4), this value for AQq is 
given by 


AQ 


0 ‘ 



(6.4.1) 


To speed up the iterative process involved in the global solution of equation (3.3.4), the 
local incremental changes given by equation (6.4.1) are multiplied by the input value of 
the overrelaxation factor O RELAX. Thus, 

A Q 0 = ^jCORELAX) (6.4.2) 


where 


1.0 < ORELAX < 2.0 (6.4.3) 

The optimum value of ORELAX for the shortest running time probably varies somewhat 
with the boundary conditions of the problem. A preliminary investigation indicated an 
approximate value of 1.35; however, the user is encouraged to try other input values. 

6.5 Accuracy (EPS and EPSR) . - The input values of EPS and EPSR determine 
the accuracy of various iterative processes in the program. The input value of EPSR 
relates to the solution of the governing equation (3.3.4) by finite-difference methods. It 
is the maximum allowable value of the residual error St at any point in the flow field 
after the iterative solution has been completed globally for a given (fixed) set of the 
coefficients Cq, Cq,..., C 5 . The input value of EPS relates to other iterative processes. 

The program uses double precision. The input values of EPS and EPSR used for 
the examples in this report were both 0.000005 and occasionally an order of magnitude 
less. Because of the dimensionless form under which the solutions are obtained, the 
magnitudes of EPS and EPSR are independent of the size of the flow field. However, 
for comparable accuracy in solving the governing equation (3.3.4), the greater the number 
of grid points on a potential surface, the smaller the input value of EPSR, because the 
smaller will be the dimensionless grid spacings a^,...,ag in the coefficients of equation 
(3.3.4). 



7.0 INPUT TO PROGRAM 


The two major inputs to the program are the shape of the upstream boundary with 
its associated grid and the velocity distribution on the lateral boundary of the duct. These 
inputs are discussed in detail in the next two sections, after which a formatted, line-by- 
line description of the complete input is given. This latter section constitutes a users 
guide for preparation of the input. 

7.1 Upstream boundary shape and associated grid . - The shape of the upstream 
boundary is specified by the coordinate points of its contour on the y,z plane for x = 0. 
These coordinate points are located at every intersection of the contour with a specified 
(input) grid of YG(J) and ZG(K) lines in the Y,Z plane. Thus, 



The coordinate points [YC(1), ZC(l)],..., (YC(NCP), ZC(NCP)] are numbered counter - 
clockwise and consecutively from 1 to NCP, where the maximum allowable value of 
NCP is 200. The starting point is arbitrary, except for cases using planar symmetry 
(ISYM equal to 2 or 3), which are discussed later. The YG(J) and ZG(K) grid lines, at 
which the contour coordinate points occur, are numbered 1 to NJ and 1 to NK, 
respectively, where the maximum allowable values of NJ and NK are 21 and 36, 
respectively. 

The spacing (a2, and ag) of the grid lines is arbitrary, except that at least 

three internal grid points (i.e., intersections of grid lines) must lie along every internal 
grid line segment bounded by the contour and at least two external grid points must lie 
along every external grid line segment bounded by the contour. It is also prudent to keep 
the grid spacings as nearly constant as the contour shape and other considerations permit 
and not too different from the a^ and a. 4 , spacing in the cp direction (second fig. in 
section 3.2). Of course, grid size affects running time. Doubling the grid spacing on 
potential surfaces, but leaving the spacing between potential surfaces unchanged, 
decreased CPU time by more than 75 percent in the following examples. 
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Fine grid; ITER ■ 10; 370/3033 CPU time, 132.01 min 


Coarse grid (spacing doubled); ITER ■ 10; 370/3033 CPU time, 29. 64 min 


The shape of the contour is completely arbitrary (but see section 4.1) except that 
(see following figure) (1) every contour point must be the end point of at least one internal 
grid line, (2) any interior straight line drawn between two contour points that are not 
adjacent must cut at least one grid line, and (3) unless a contour point lies on a grid point, 
it must be at least 10 times the input value of EPS away from any interior grid point. 
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Because of the approximate nature of the solutions (section 4.0), a "kink" may 
develop in the duct boundary if an external grid point is too close to two adjacent contour 
points along the boundary, particularly if the upstream boundary configuration is too 
complex or convoluted. This "kink" occurs because the x,y,z coordinates at boundary 



point A are obtained by extrapolating from the corresponding values at the interior 
points 1, 2, and 3, whereas the coordinates at boundary point B are obtained by extrap- 
olating from the interior points a, b, and c. The "kink" is most easily eliminated by 
shifting one of the two grid lines so that the points A and B come together as shown in 
the following figure. 
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For planar symmetry solutions (section 6.2), only one of the two symmetrical halves 
of the upstream boundary configuration is used. If the plane of symmetry is normal to the 
Y axis (ISYM = 2) as shown, the left half must be used, and the contour points are num- 
bered in the counterclockwise direction from 1 to NCP, starting at the plane of 
symmetry. 



If the plane of symmetry is normal to the Z axis (ISYM = 3), the lower half must be 
used, and the contour points are numbered in the counterclockwise direction from 1 to 
NCP starting at the plane of symmetry as shown. Note that for ISYM equal to both 2 



and 3, the input location ( JX,KX) of the primary streamline must be on the plane of 
symmetry. 

These various upstream boundary contours in physical X,Y,Z space are also the 
shapes of all potential surfaces in transformed <p,t|r,Ti space, because paired values of Y 
(equals fy) and Z (equals p) are constant along streamlines (section 3.2). 

7.2 Prescribed velocity distribution on surface of duct . - The velocity distribution 
on the lateral boundary of the duct could be specified in a perfectly general, continuous 
way (but see section 4.1) at each of the NCP coordinate points along the boundary 
contour (section 7.1) for each of the NI potential surfaces from the upstream boundary 
(I = 1) to the downstream boundary (I = NI), where the maximum allowable value for both 
NCP and NI is 200. Because this is a large amount of input data (200 x 200), for conven- 
ience, in program DIN3D1, the velocity distribution on the lateral surface is specified by 
two components. First, the distribution of velocity QP(I) is specified as a function of 
distance SP(I) along the principal streamline ( input value of IP; figs, in section 7.1). 

Thus, 
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where QP(I) is expressed as a ratio of the upstream velocity, and the dimensional unit 
for SP(I) is the same as for YG(J) and ZG(K) (section 7.1). The velocity QP(I) is 
constant in the upstream and downstream regions, which regions should normally be at 
least two hydraulic diameters of their respective flow areas in extent. These regions of 
constant velocity on the lateral boundary are required to justify the assumption of con- 
stant velocity over the upstream and downstream flow areas. 

Second, but only if the input value of option IVEL is 1, the velocity variation DQ 
around the contour of each potential surface, which contour in cp,i|x,r| space is the same 
as the upstream contour (section 7.1), is specified by 



In this figure, DQAMP(I) is the amplitude (plus or minus) of the velocity variation DQ 
and P is the decimal fraction of the distance around the contour. The XP(I) value of P 
locates the principal streamline IP relative to the velocity variation with P. The vari- 
ation in velocity with P is, therefore, specified for all potential surfaces by DPOKI), 

DP 12(1), DP23(I), DQAMP(I), and XP(I) as functions of SP(I) from I = 1,...,NI. As for the 
distribution of QP(I) in the previous paragraph, the distributions of these parameters 
should also be constant in the upstream and downstream regions, and the values of 
DQAMP(I) must be zero. 

In the regions of P defined by DPOKI) and DP23(I) in the figure, the velocity 
variation DQ is given by the cubic equation 

DQ = a + bP + cP 2 + dP 3 (7.2.1) 
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where the four coefficients a, b, c, and d are fixed by the four conditions 

d(DQ) 

=0 at the two end points 

DQ = 0 at one end point 

DQ = DQAMP(I) at the other end point 

To simplify the input further, values of these parameters, as well as of QP(I), need 
not be specified at all values of I, but only at NP values, where 

NI = (NP - 1) NSD + 1 (7.2.2) 

in which NSD is the specified (input) number of equal subdivisions between adjacent, 
specified NP values of the parameters. The program assumes linear variations in the 
parameters between the specified values. 

For the "equilibrium" velocity distributions described in appendix A ( input values of 
option IVEL equal to 2 or 3), in addition to the prescribed velocity QP(I) as a function 
of distance SP(I) along the primary streamline, only the amplitude DQAMP(I) of the 
velocity variation DQ (see previous fig.) is specified. The parameters DPOl(I), DP 12(1), 
DP23(I), and XP(I) must be omitted. Also, if the input value of ISYM is 2, the input 
value of IVEL must not be 3; and if the input value of ISYM is 3, the input value of 
IVEL must not be 2. 

Option IVEL = 4 can be used only with option ISYM equal to 2 or 3 (planar symme- 
try cases). Here, only the parameters DP23(I) and XP(I) must be omitted, it being 
understood that for planar symmetry 

DP23(I) = DPOl(I) (7.2.3) 


and 


XP(I) = 0.5 + DPOl(I) + 0.5 x DP 12(1) (7.2.4) 

As for optional input IFLUID (section 6.1), provision is also made in subroutine 
PUTIN for adding new types of option IVEL by additions to the code, and of course, 
appropriate additions to the code must also be made in subroutine VELD. 

7.3 Line-by-line input for program DIN3D1 . - This section should be used when 
preparing the formatted, line-by-line input for program DIN3D1. It is also recommended 
that sections 7.1 and 7.2 be reviewed before starting. 


Line 1 - FORMAT(2QA4) 

TITLE title (center on field of 80 characters) 

Line 2 - FQRMAT(20A4) 

SUBT1 first subtitle (center on field of 80 characters) 

Line 3 - FORMAT(2QA4) 

SUBT2 second subtitle (center on field of 80 characters) 
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Line 4 - F0RMAT(7I1Q) 

IFLUID option equals 1 for incompressible flow, 2 for perfect gas (section 6.1) 

ISYM option equals 1 for complete flow field, 2 for half flow field with planar 

symmetry about y plane, and 3 for half flow field with planar symmetry 
about z plane (sections 6.2 and 7.1) 

IGRID option equals 1 for Cartesian YG(J),ZG(K) grid at upstream boundary 

(only option) (section 7.1) 

IVEL option equals 1 for standard, two -component, parametric method of 

specifying velocity distribution on lateral boundary of duct (section 7.2), 

2 for "equilibrium" velocity distribution with turn in y plane (sections 
4.1 and 7.2 and appendix A), 3 for "equilibrium" velocity distribution with 
turn in z plane (sections 4.1 and 7.2 and appendix A), and 4 for cases 
with ISYM values of 2 or 3 only (section 7.2) 

IPLOT option equals zero for no graphics output, 1 for three-dimensional 

graphics output (section 6.3) 

ISOLV option equals zero if built-in criteria for successful solution are not 

used, 1 if criteria are used (section 4.4) 

ISPACE option equals zero if grid spacings (second fig. in section 3.2) for all 

internal grid points are not printed in output (space-saving option), 1 if 
spacings are printed 

Line 5 - FORMAT(SIIO) 

ITERMX maximum number of major iterations (ITER) allowed; value depends on 

circumstances of case involved (sections 3.3, 4.4, 4.5, and 5.1) 

ITMAX maximum number of IT iterations allowed (each iteration involves 

entire flow field; values of coefficients Cq, Cq> C ^,...,Cg in equation 
(3.3.4) are unchanged for all IT iterations; recommended value is 250) 
(sections 3.3 and 5.1) 

ITXMAX maximum number of passes allowed for iterative, finite-difference 

solution of governing differential equation (3.3.4) on a given potential 
surface <p; recommended value is 100 (section 5.19) 

ICONX maximum number of ICON iterations, in main program, on coefficients 

of governing equation (3.3.4); recommended value is 4 and cannot be less 
than 3 (section 5.1.1) 

NTRY number of iterations, in subroutine VARI, on values of direction cosines; 

recommended value is 3 and cannot be less than 2 

Line 6 - FQRMAT(6F10.3) 

CAVD coefficient for averaging new values of derivatives of direction cosines 

with previous values; recommended value is 0.5 (section 5.2.2) 
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CAVN coefficient for averaging new values of continuity parameters A and B 

with previous values; recommended value is 0.5 or less (section 5.2.2 and 
appendix C) 

CAVP coefficient for averaging values (obtained by two methods) of x,y,z 

coordinates at each internal grid point; recommended value is 0.5 or less 
(sections 3.6, 4.2(constraint 1), and 4.6) 

CAVX coefficient for averaging new values of direction cosines with previous 

values; recommended value is 0.2 (section 5.2.2) 

CAVY coefficient for averaging new values of coefficient C£, in the governing 

equation (3.3.4), with previous values; recommended value is 0.5 or less 
(section 5.2.2) 

CAVZ coefficient for averaging new values of THET (cosine of "distortion" 

angle 9) with previous values; recommended value is 0.5 or less (section 
5.2.2) 

[If IFLUID = 1, incompressible flow, go to line 8.] 

Line 7 - FQRMAT(2F10.4) 

AMU upstream Mach number (section 6.1) 

GAM ratio of specific heats (section 6.1) 

Line 8 - FORMAT(2UO) 

JX value of J for primary streamline (fig. in section 3.5; sections 3.6, 5.2, 

and 7.1) 

KX value of K for primary streamline (fig. in section 3.5; sections 3.6, 5.2, 

and 7.1) (For input values of ISYM equal to 2 and 3, the primary 
streamline (JX,KX) must lie on the plane of symmetry.) 

Line 9 - FQRMAT(3I10) 

NJ number of YG(J) grid lines; maximum value is 21 (section 7.1) 

NK number of ZG(K) grid lines; maximum value is 36 (section 7.1) 

NCP number of contour coordinate points around upstream potential surface; 

maximum value is 200 (sections 6.3 and 7.1) 

Line 10 - FQRMAT(8F10.6) 

YG(J) NJ values of Y grid lines; same dimensional unit of length used for SP 

on line 15 ; 6-decimal accuracy recommended (sections 3.2, 3.5, and 7.1) 

Line 11 - FQRMAT(8F10.6) 

ZG(K) NK values of Z grid lines; same dimensional unit of length used for SP 

on line 15 ; 6-decimal accuracy recommended (sections 3.2, 3.5, and 7.1) 
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Line 12 - FQRMAT(8F10.6) 

YC(IX) NCP values of Y for coordinate points along boundary contour starting 

at contour point 1 (which has arbitrary location for ISYM = 1 but must 
lie on the plane of symmetry for ISYM equal to 2 or 3); contour points 
must be read sequentially in counterclockwise direction; same 
dimensional unit of length used for SP on line 15 ; 6-decimal accuracy 
recommended (sections 3.2, 3.5, and 7.1) 

Line 13 - FQRMAT(8F10,6) 

ZC(IX) NCP values of Z for coordinate points along boundary contour starting 

at contour point 1; see YC(IX), above, for further comments 

Line 14 - FORMAT(3I10) 

IP contour coordinate point corresponding to principal streamline; for 

ISYM equal to 2 or 3, IP must equal 1; for "equilibrium" velocity 
distributions (IVEL equal to 2 or 3), see appendix A (first and last figs, in 
section 7.1) 

NP number of stations (velocity potential surfaces) at which parameters are 

specified for velocity distribution on lateral surface; quantity 
(NP - 1)NSD + 1 must not exceed 200 (section 7.2) 

NSD number of subdivisions between each of the above NP stations 

(section 7.2) 

Line 15 - (NP lines, one for each station), FQRMAT(6F10.5, F10.6) 

QP(I) velocity (ratio) distribution along principal streamline (IP), expressed as 

ratio of upstream velocity; 5-decimal accuracy recommended 
(section 7.2) 

SPG) distance along principal streamline; any unit of length permitted, and 

value of SPG) at upstream boundary need not be 0.0 (section 7.2) 

DPOl(I) percent of contour length (second fig. in section 7.2); omit if IVEL 

equals 2 or 3 (section 7.2) 

DP12(I) percent of contour length (second fig. in section 7.2); omit if IVEL 

equals 2 or 3 (section 7.2) 

DP23(I) percent of contour length (second fig. in section 7.2); omit if IVEL 

equals 2, 3, or 4 (section 7.2) 

DQAMP(I) amplitude (second fig. in section 7.2) of velocity variation DQ around 
contour of potential surface; velocity expressed as ratio of upstream 
velocity; DQAMP(I) may be positive or negative; DQAMP(I) must be 0.0 
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in upstream and downstream regions of duct (first fig. in section 7.2; 
section 7.2) 


XP(I) 

Line 16 - 
IA 

IZ 

NPS 

NSL 

Line 17 - 
IPS(I) 

Line 18 - 
ISL(I) 

Line 19 - 
EPS 

EPSR 

ORELAX 


location (percent of contour length) of principal streamline relative to 
velocity variation around contour of potential surface (second fig. in 
section 7.2); omit if IVEL equals 2, 3, or 4 (section 7.2) 

FQRMAT(4I10) 


I value of initial potential surface for which output data are printed in 
table I 

I value of final potential surface for which output data are printed in 
table I; IZ £NI (total number of potential surfaces; section 6.3); if 
IZ < IA, table I is omitted in printout 

number of potential surfaces for which output data at boundary contour 
points are printed in table II (and saved for three-dimensional graphics if 
input value of IPLOT is 1); NPS < 200; if NPS = 0, table II is omitted 
and input line 17 is skipped (section 6.3) 

number of boundary-surface streamlines for which output data are 
printed in table III (and saved for three-dimensional graphics if IPLOT 
is 1); NSL < 200; if NSL = 0, table III is omitted and input line 18 is 
skipped (section 6.3) 

FQRMAT(8I10) 


NPS values of the I values of potential surfaces for which output data 
are printed in table II; numbered sequentially, starting from lowest 
value, but numbers can be skipped (section 6.3) 

FQRMAT(8I10) 


NSL values of the I values of boundary contour points for which 
streamline data are printed in table III; numbered sequentially, starting 
from lowest value, but numbers can be skipped (section 6.3) 

FQRMAT(2F10.7,F10.4) 

standard maximum allowable error in various iterative procedures; 
recommended value is 0.000005 (section 6.5) 

maximum allowable value of residual error 0t in finite-difference 
solutions of equation (3.3.4); recommended value is 0.000005 (sections 
3.3, 5.19, and 6.5) 

overrelaxation factor (sections 5.19 and 6.4) 
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7.4 Sample printout of input data . - Program DIN3D1 prints out the input data in 
the same order in which they are read in. A sample printout of the input data follows. 


PROGRAM DIN3D1 

DESIGN OF THREE-DIMENSIONAL INTERNAL FLOW FIELDS 
FOR ARBITRARY PRESCRIBED VELOCITY DISTRIBUTIONS 
ON LATERAL BOUNDARY SURFACE 


CASE NO. V 
ELBOW C2 

QD/QU =2.0 M-UP =0.4 


INPUT DATA 
OPTIONS 

I FLUID ISYM IGRID IVEL IPLOT ISOLV ISPACE 

2 112 111 

INPUT DATA FOR LIMITS ON VARIOUS ITERATION CYCLES 

MAX ITER MAX IT MAX ITX ICONX NTRY 

ITERATIONS ITERATIONS ITERATIONS ITERATIONS ITERATIONS 

12 250 100 4 3 

INPUT DATA FOR VARIOUS DAMPING COEFFICIENTS 

CAVD CAVN CAVP CAVX CAVY CAVZ 

0.500 0.500 0.200 0.200 0.500 0.500 

INPUT DATA FOR PERFECT GAS CIFLUID = 2) 

UPSTREAM RATIO OF 
MACH NO. SPEC HTS 

0 . A0 00 1.4Q00 

INPUT DATA FOR PRIMARY STREAMLINE 

J-VALUE K-VALUE 
(JX) CKX) 
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INPUT 

HO, 

GRID 


INPUT 


INPUT 


INPUT 


DATA FOR GRID SYSTEM 

OF Y NO. OF Z 

LINES GRID LINES 

13 II 


OH UPSTREAM BOUNDARY 

NO. OF POINTS 
ON BOUNDARY 

28 


SURFACE (IGRI'D = I) 


DATA FOR UPSTREAM GRID (CONTINUED) 


J 

Y-VALUE 
OF GRID 

1 

1.000000 

2 

1 .290635 

3 

2.313019 

4 

3.909010 

5 

4.704505 

6 

5.500000 

7 

6.295495 

8 

7.090990 

9 

7 . 886485 

10 

8.681930 

11 

10.272971 

12 

11.300355 

13 

11.590990 


DATA FOR UPSTREAM GRID (CONTINUED) 

Z-VALUE 
K OF GRID 


1 1.000000 

2 1.290635 

3 2.31S0I9 

4 3 . 9 0 9 0 1 G 

5 4.704505 

6 5.500QGG 

7 6.295495 

3 7.090990 

9 8.681981 

10 9.7 09365 

11 10.000000 


DATA FOR UPSTREAM GRID (CONTINUED) 

Y-VALUE Z-VALUE 

I OF CONTOUR OF CONTOUR 


1 

5 . 50 0 O 0 G 

1.000000 

2 

6.295495 

I. 000000 

3 

7.090990 

1.000000 

4 

7.886485 

1.070871 

5 

8.631930 

1.290635 

6 

10.272971 

2.318019 

7 

11.300355 

3.909010 

8 

11.520119 

4.704505 

9 

11.590990 

5 . 5 0 0 G 0 0 

10 

11.520119 

6.295495 

11 

11 . 30 0355 

7.090990 

12 

10.272971 

8. 681931 

13 

8.681980 

9.709365 

14 

7.886435 

9.929129 

15 

7.090990 

10.000000 

16 

6.295495 

10.000000 

17 

5.500000 

10.000000 

13 

4.704505 

9.929129 

19 

3.909010 

9.709365 
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20 

2.318019 

8.681931 

21 

1.290635 

7.090990 

22 

1.070871 

6.295495 

23 

1.000000 

5.500000 

24 

1.070371 

4.704505 

25 

1.290635 

3. 909010 

26 

2.318019 

2.318019 

27 

3.9Q9G1G 

1.290635 

28 

4.704505 

1.070871 


INPUT DATA FOR VELOCITY DISTRIBUTION ON LATERAL BOUNDARY SURFACE CIVEL = 1) 



PRINCIPAL NO. 

GF SPEC. 

NO. OF 

STREAMLINE 

STATION'S SUBDIVISIONS 


16 

29 

1 


VEL. ON 

DISTANCE 



PRINCIPAL 

ALONG 

DEL-Q 

I 

STREAMLINE 

STREAMLINE 

AMPLITUDE 

1 

1.000000 

0 . 0 0 COG 0 

0.000000 

o 

c. 

1.000000 

1.500000 

0.000000 

3 

1.000000 

3.000000 

0.000000 

4 

1.000000 

4.500000 

0.000000 

5 

1 . G 0 0 0 0 0 

6 . C 0 0 0 0 0 

0 . 000000 

6 

l.OOOGOO 

7.500000 

O.QOOOOQ 

7 

1.000000 

9.000000 

0.000000 

8 

1.019700 

10.500000 

“*0.019700 

9 

1.074100 

12.000000 

“0 . 074100 

10 

1.156300 

13.500000 

-0.156300 

11 

1.259300 

15. 00G00Q 

-0 .259300 

12 

1.376200 

16.500G00 

-0.356500 

13 

1.500000 

18. D 0 0 0 D G 

-0.425900 

14 

1.623800 

19.500000 

-0.467500 

15 

1.740700 

21.000000 

-0.481400 

16 

1.843800 

22.500000 

-0.467600 

17 

1.925900 

24. 000000 

-0.425900 

13 

1.9S0300 

25.500000 

“0 . 356 50 0 

19 

2.000000 

27.000000 

-0.259300 

20 

2.000000 

28.500000 

-0.156200 

21 

2.000000 

30.000000 

-0.074100 

22 

2.000000 

31.500000 

-0.019700 

23 

2.000000 

33.000000 

0.000000 

24 

2.000000 

34.500000 

0.000000 

25 

2 . GO 0 0 0 0 

36 .000000 

0.000000 

26 

2.000000 

37.500000 

0 .000000 

27 

2.000000 

39.000000 

0. 000000 

23 

2 . 0 Q 0 0 0 0 

40.500000 

0.000000 

29 

2.000000 

42.000000 

0 .000000 


INPUT DATA FOR PRINT FORMAT 


MIN I-VALUE 
OF POT SURF 
(TABLE I) 


MAX I ~ VALUE 
OF POT SURF 
(TABLE I) 


NO OF POT 
SURFACES 
(TABLE II) 


NO OF 
STREAMLINES 
(TABLE III) 


1 


29 


29 


28 
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INPUT DATA FOR POTENTIAL SURFACES (TABLE II) 


NUMBER (NFS) OF POTENTIAL SURFACES = 29 
I-VALUES OF POTENTIAL SURFACES: 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

IS 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 



INPUT DATA FOR STREAMLINES (TABLE III) 

NUMBER (NSL) OF STREAMLINES = 28 

I-VALUES OF CONTOUR POINTS THROUGH WHICH STREAMLINES PASS: 


1 

2 

3 

4 

5 

6 

7 

3 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 




INPUT DATA RELATED TO ACCURACY OF CALCULATIONS 


EPS 

0.0000005 


EPS-R 

0.0000005 


G-RELAX 

1.3500 
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8.0 OUTPUT FROM PROGRAM 


The five major outputs from the program are (1) an intermediate printout generated 
as the solution progresses; (2) output table I, with data at the internal grid points for the 
selected range (IA to IZ) of potential surfaces in x,y,z space; (3) output table II, with 
data around the contours of selected potential surfaces in x,y,z space; (4) output 
table III, with data along selected streamlines over the full range of potential surfaces 
(I equals 1 to NI); and (5) output data to tape or disk for three-dimensional graphics, 
provided that the input value of IPLOT is 1. Before printing these five outputs, the 
program prints out (1) the maximum Mach number along the principal streamline and 
(2) provided that the input value of ISPACE is 1, the values of the six grid spacings 
ai,...,a$ (second fig. in section 3.2) for all of the internal grid points. 

8.1 Intermediate printout . - For each pass IT through the entire flow field for 
every major iteration ITER (with unchanged values of the coefficients Cq, Cq» C]_,...,C$ 
in the governing equation (3.3.4)), the intermediate printout gives the magnitude and 
location (I-MAX, J-MAX, and K-MAX) of the absolute value of the maximum residual 
error 31 encountered during the pass. For a given value of ITER, after convergence (Si < 
EPSR) or after IT becomes greater than the input value of ITMAX, the program prints 
the exit flow area EXFLAR (section 5.12) and its error ERRAR (section 5.1), expressed 
as a decimal fraction, for each of the ICON iterations (section 5.1.1). Also, any 
intermediate messages regarding, for example, the counters ICX, ICY, and ITH 
(section 5.2.1) are printed. A sample page of intermediate printout follows. 


INTERMEDIATE PRINTOUT 


ITERATION 

ITERATION 

MAX RES 




ITERATION 

EXIT FLOW 

CORRECT FLOW 


ITER 

IT 

IN IT 

I-MAX 

J-MAX 

K-MAX 

ICON 

AREA (DIM) 

AREA (DIM) 

ERROR 

1 

1 

0.8707286 

8 

7 

10 





I 

2 

0.3962490 

8 

7 

8 

... 

... 



1 

3 

U. 2352487 

7 

7 

7 


... 



1 

4 

0.1334014 

7 

7 

7 





1 

5 

0.0794694 

6 

7 

7 


... 



1 

6 

0.0498247 

6 

7 

6 


. . . 



1 

7 

0.0313226 

6 

7 

6 





1 

8 

0.0192473 

6 

7 

6 





1 

9 

0.0120385 

5 

7 

6 





1 

10 

0.0073951 

5 

7 

6 





1 

11 

0.0044359 

5 

7 

6 





1 

12 

0.0027005 

4 

7 

6 





1 

13 

0.0016138 

4 

7 

6 





1 

14 

0.0009475 

4 

7 

6 





1 

15 

0.0005484 

4 

7 

6 





1 

16 

0.0003137 

4 

7 

6 





1 

17 

0.0001780 

4 

7 

6 





1 

18 

0.0001002 

4 

7 

6 





1 

19 

0.0000562 

4 

7 

6 





1 

20 

0.0000315 

4 

7 

6 





1 

21 

0.0000176 

4 

7 

6 





1 

22 

0.0000098 

4 

7 

6 





1 

23 

0.0000055 

4 

7 

6 





1 

24 

0.0000031 

4 

7 

6 





1 

25 

0.0000017 

3 

7 

6 





1 

26 

0.0000010 

4 

7 

6 





1 

27 

0.0000005 

7 

7 

6 

* 




I 

27 

. . . 

... 

. . . 

. . . 


49.3908 

49.3932 

-0.0000 

1 

27 

. . . 

... 


. . . 

2 

49.3927 

49.3932 

-0 . 000 0 

1 

27 

. . . 

... 

. . . 

... 

3 

49.3935 

49.3932 

0.0000 

1 

27 





4 

49.3939 

49.3932 

0.0000 

2 

1 

0.0682211 

12 

7 

2 





2 

2 

0.0395768 

12 

7 

6 





2 

3 

0.0248684 

12 

7 

6 





2 

4 

0.0154618 

11 

7 

6 





2 

5 

0.0096363 

10 

7 

6 





2 

6 

0.0060232 

10 

7 

6 





2 

7 

0.0039346 

9 

7 

6 





2 

8 

0.0024947 

8 

7 

6 





2 

9 

0.0016270 

8 

7 

6 



. . . 
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8.2 Output table I . - Output table I gives the values of variables at all internal grid 
points for potential surfaces in the (input) range from IA to IZ. The headings in output 
table I are as follows: 

I, J,K grid-point indices in directions of increasing velocity potential PHI, 

stream function PSI, and stream function ETA, respectively 

PHI, PSI, at grid point (I,J,K), values of velocity potential and two stream 

ETA functions, respectively (sections 3.0 to 3.2) 

X(DIM),Y(DIM), at grid point (I, J,K), values of X,Y,Z coordinates, expressed in same 
Z(DIM) dimensional unit as input values of SP(I) (sections 3.5 and 3.6) 

Q/Q-UP at grid point (I, J,K), value of local velocity divided by upstream velocity 

(section 3.3) 

MACH NO. at grid point (I, J,K), value of local Mach number 

RO/RO-UP at grid point (I, J , K), value of local static density divided by upstream 

static density (section 3.3) 

P/P-UP at grid point (I, J,K), value of local static pressure divided by upstream 

static pressure. For incompressible flow (IFLUID = 1), P/P-UP is 
defined as local difference between total and static pressure divided by 
the same difference at upstream boundary (which definition is equivalent 
to square of Q/Q-UP) (section 3.3)) 

SIN(THET) sine of "distortion" angle 0 (sections 3.1, 3.3, and 4.2 (constraint 5)) 

COS(ALl),..., at grid point (I,J,K), values of direction cosines of three unit vectors 

COS(GM3) ej_, e2> and e% (sections 3.1, 3.4, and 4.2) 

A,B at grid point (I,J,K), values of continuity parameters (eqs. (3.3.2) and 

(3.3.3) and sections 3.4 and 5.2) 
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A sample page of output table I resulting from the sample input in section 7.4 follows 


OUTPUT TABLE NO. I (INTERNAL GRID POINTS) 

I J K PHI PSI ETA X( DIM) Y(DIM) Z(DIM) Q/Q-UP MACH NO RO/RO-UP P/P-UP 5INCTHET) 

C0SCAL1) C0SCBT1) C0SCGM1) C0SCAL2) COS ( BT2 ) C0SCGM2) C0SCAL3) C0SCBT3) C0SCGM3) A B 

14 10 9 2.51638 0.87683 0.87683 19.7331 7.4296 8.1175 1.5330 0.6269 0.8687 0.8212 1.0000 

0.9757 -0.0309 0.2171 0.0316 0.9995 0.0004 -0.2170 0.0065 0.9762 0.3018 0.7285 

14 5 10 2.51638 0.42284 0.99409 19.5401 3.8690 8.9232 1.6026 0.6578 0.8509 0.7977 1.0000 

0.9783 0.0205 0.2060 -0.0209 0.9998 -0.0003 -0.2059 -0.0045 0.9786 0.7891 0.6952 

14 6 10 2.51638 0.51363 0.99409 19.5290 4.5824 8.9232 1.6026 0.6577 0.8510 0.7978 1.0000 

0.9785 0.0103 0.2060 -0.0105 0.9999 -0.0001 -0.2060 -0.0023 0.9786 0.7890 0.6954 

14 7 10 2.51638 0.60^43 0.99409 19.5252 5.2955 8.9233 1.6026 0.6577 0.8510 0.7978 1.0000 

0.9786 0.0000 0.2060 -0.0000 1.0000 -0.0000 -0.2060 -0.0000 0.9786 0.7889 0.6955 

14 8 10 2.51638 0.69523 0.99409 19.5290 6.0086 8.9232 1.6026 0.6577 0.8510 0.7978 1.0000 

0.9785 -0.0103 0.2060 0.0105 0.9*?9 0.0001 -0.2060 0.0023 0.9786 0.7890 0.6954 

14 9 10 2.51638 0.78603 0.99409 19.5401 6.7220 8.9232 1.6026 0.6578 0.8509 0.7977 1.0000 

0.9783 -0.0205 0.2060 0.0209 0.9998 0.0003 -0.2059 0.0045 0.9786 0.7891 0.6952 


15 5 2 2.80440 0.42284 0.03317 23.4195 3.9438 2.3670 1.2708 0.5134 0.9277 0.9003 1.0000 

0.9258 0.0180 0.3775 -0.0196 0.9998 0.0003 -0.3775 -0.0063 0.9260 0.8083 0.9264 

15 6 2 2.80440 0.51363 0.03317 23.4093 4.6195 2.3671 1.2709 0.5134 0.9277 0.9003 1.0000 

0.9259 0.0090 0.3776 -0.0097 1.0000 -0.0000 -0.3776 -0.0032 0.9260 0.3083 0.9263 

15 7 2 2.80440 0.60443 0.03317 23.4060 5.2955 2.3672 1.2709 0.5134 0.9277 0.9003 1.0000 

0.9260 0.0000 0.3776 -0.0000 1.0000 -0.0000 -0.3776 -0.0000 0.9260 0.3083 0.9262 

15 8 2 2.80440 0.69523 0.03317 23.4093 5.971+ 2.3671 1.2709 0.5134 0.9277 0.9003 1.0000 

0.9259 -0.0^90 0.3776 0.0097 1.0000 0.0000 -0.3776 0.0032 0.9260 0.3083 0.9263 

15 9 2 2.80440 0.78603 0.03317 23.4195 6.6472 2.3670 1.2708 0.5134 0.9277 0.9003 1.0000 

0.9258 -0.0180 0.3775 0.0196 0.9998 -0.0003 -0.3775 0.0063 0.9260 0.8083 0.9264 

15 4 3 2.80440 0.33204 0.15044 23.0727 3.2591 3.2840 1.3129 0.5313 0.9190 0.8885 1.0000 

0.9318 0.0275 0.3620 -0.0297 0.9996 0.0006 -0.3619 -0.0088 0.9322 0.5046 0.8928 

15 5 3 2.80440 0.42284 0.15044 23.0555 3.9377 3.2842 1.3130 0.5314 0.9192 0.8887 1.0000 

0.9319 0.0182 0.3621 -0.0197 0.9998 0.0004 -0.3621 -0.0059 0.9321 0.8043 0.8930 

15 6 3 2.80440 0.51363 0.15044 23.0453 4.6164 3.2843 1.3131 0.5314 0.9192 0.8888 1.0000 

0.9321 0.0091 0.3622 -0.0098 1.0000 0.0001 -0.3622 -0.0030 0.9321 0.3044 0.8929 

1573 2.80440 0.60443 0.15044 23.0419 5.2955 3.2843 1.3131 0.5314 0.9192 0.8888 1.0000 

0.9321 0.0000 0.3622 -0.0000 1.0000 -0.0000 -0.3622 -0.0000 0.9321 0.3044 0.8928 

15 8 3 2.80440 0.69523 0.15044 23.0453 5.9745 3.2843 1.3131 0.5314 0.9192 0.8888 1.0000 

0.9321 -0.0091 0.3622 0.0098 1.0000 -0.0001 -0.3622 0.0030 0.9321 0.3044 0.8929 

15 9 3 2.80440 0.78603 0.15044 23.0555 6.6533 3.2842 1.3130 0.5314 0.9192 0.8887 1.0000 

0.9319 -0.0182 0.3621 0.0197 0.9998 -0.0004 -0.3621 0.0059 0.9321 0.3043 0.8930 

15 10 3 2.80440 0.87683 0.15044 23.0727 7.3319 3.2840 1.3129 0.5313 0.9190 0.8885 1.0000 

0.9318 -0.0275 0.3620 0.0297 0.9996 -0.0006 -0.3619 0.0088 0.9322 0.3046 0.8928 

15 3 4 2.80440 0.15044 0.33204 22.6111 1.8814 4.6653 1.3828 0.5614 0.9037 0.8678 1.0000 

0.9389 0.0473 0.3408 -0.0505 0.95S7 0.0004 -0.3405 -0.0143 0.9401 0.7969 0.8427 

15 4 4 2.80440 0.33204 0.33204 22.5550 3.2461 4.6653 1.3828 0.5614 0.9042 0.8685 1.0000 

0.9396 0.0279 0.3410 -0.0298 0.9996 0.0003 -0.3409 -0.0085 0.9401 0.7971 0.8425 

15 5 4 2.80440 0.42284 0.33204 22.5376 3.9290 4.6653 1.3828 0.5614 0.9044 0.8687 1.0000 

0.9399 0.0185 0.3410 -0.0197 0.9998 0.0002 -0.3410 -0.0056 0.9400 0.7972 0.8424 

15 6 4 2.80440 0.51363 0.33204 22.5273 4.6122 4.6654 1.3829 0.5614 0.9044 0.8688 1.0000 

0.9400 0.0092 0.3411 -0.0098 1.0000 0.0001 -0.3411 -0-0028 0.9400 0.7973 0,8423 
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8.3 Output table n . - Output table II gives the values of variables along the 
boundary contour of potential surfaces selected by the input. The headings in output table 
II are as follows: 


I index number for potential surface (constant PHI) 

ICP index number for contour point along boundary of potential surface I 

X-CP(DIM), at contour point ICP of potential surface I, values of X,Y,Z 
Y-CP(DIM), coordinates, expressed in same dimensional units as input values 
Z-CP(DIM) of SP(I) (sections 5.6 and 6.3) 


Remaining headings for table II are defined under output table I (section 8.2). A sample 
page of output table n resulting from the sample input in section 7.4 follows. 


OUTPUT TABLE NO. II 

(COORDINATE 

POINTS ALONG 

CONTOURS 

OF SELECTED 

POTENTIAL 

surfaces: 

I 

ICP 

X-CP(DIM) 

Y-CP(DIM) 

Z-CP ( DIM) 

Q/Q-UP 

MACH NO 

RO/R0-UP 

P/P-UP 

17 

16 

23.7339 

5.2955 

10.6006 

1.9259 

0.8061 

0.7972 

0.7281 

17 

17 

23.7365 

4.6218 

10.6002 

1.9259 

0 .8061 

0.7972 

0,7281 

17 

18 

23.7679 

3.9486 

10.5500 

1.9216 

0.8041 

0.7983 

0.7295 

17 

19 

23.8549 

3.2775 

10.3952 

1.9084 

0.7978 

0.8019 

0 .7341 

17 

20 

24.2432 

1 . 9506 

9.6664 

1.8490 

0.7700 

0.8176 

0.7543 

17 

21 

24.8302 

1.1218 

8.5152 

1.7640 

0.7307 

0,8395 

0.7828 

17 

22 

25.1230 

0 . 9565 

7.9289 

1.7244 

0.7126 

0.8495 

0.7959 

17 

23 

25.4175 

0 .9150 

7.3359 

1.6865 

0.6954 

0.8589 

0.8082 

17 

24 

25.7139 

0.9903 

6.7365 

1.6502 

0.6791 

0.8678 

0.8199 

17 

25 

26 .0126 

1.1871 

6.1306 

1.6155 

0 .6635 

0.8762 

0.8310 

17 

26 

26.6218 

2.0544 

4 .9016 

1.5502 

0 . 6 345 

0.8915 

0.8515 

17 

27 

27 . 0243 

3.3597 

4.1001 

1.5108 

0.6171 

0.9005 

0.3636 

17 

28 

27.1104 

4.0064 

3.9286 

1.5026 

0.6135 

0.9024 

0.8661 

18 

1 

23.7558 

4.6607 

4.8332 

1.6238 

0.6672 

0.8742 

0.8284 

18 

2 

23.7537 

5.2955 

4.8339 

1.6238 

0.6672 

0.8742 

0 .8284 

18 

3 

28.7553 

5.9303 

4.8332 

1.6238 

0.6672 

0.8742 

0 .8284 

18 

4 

23.7292 

6.5655 

4,8844 

1.6261 

0.6683 

0.8736 

0.8276 

18 

5 

23.6372 

7.2027 

5.0450 

1.6333 

0.6715 

0.8719 

0.8254 

18 

6 

23.2043 

8.4895 

5.7998 

1.6678 

0 .6870 

0.8635 

0.8143 

18 

7 

27.5470 

9.3456 

6 . 9658 

1.7241 

0.7125 

0.8496 

0.7960 

13 

8 

27 . 2229 

9.5406 

7.5442 

1.7538 

0.7260 

0.8421 

0.7862 

18 

9 

26 . 9008 

9.6161 

8.1186 

1.7844 

0.7401 

0.8343 

0.7760 

18 

10 

26.5803 

9.5763 

8.6890 

1.8162 

0.7548 

0.8261 

0.7654 

18 

11 

26.2611 

9.4142 

9.2552 

1.8491 

0.7700 

0.8176 

0.7543 

18 

12 

25.6212 

8.5980 

10.3732 

1 . 9186 

0.8026 

0.7991 

0.7306 

18 

13 

25.1999 

7.2887 

11.0851 

1 . 9664 

0 .8253 

0.7862 

0.7141 

18 

14 

25.1064 

6.6260 

11.2367 

1.9769 

0.8303 

0.7833 

0.7104 

18 

15 

25.0736 

5.9610 

11.2864 

1 . 9803 

0.8320 

0.7824 

0.7093 

18 

16 

25.0714 

5.2955 

11.2870 

1.9803 

0 .8320 

0.7824 

0.7093 

18 

17 

25.0736 

4.6300 

11.2864 

1.9803 

0.8320 

0.7824 

0 .7093 

18 

18 

25.1064 

3.9650 

11.2367 

1.9769 

0.8303 

0.7833 

0.7104 

18 

19 

25.1999 

3.3023 

11. 0S51 

1.9664 

0.8253 

0.7862 

0.7141 

18 

20 

25.6212 

1 .9930 

10.3732 

1.9186 

0.8026 

0.7991 

0.7306 

18 

21 

26.2611 

1.1768 

9.2552 

1.8491 

0.7700 

0.8176 

0.7543 

18 

22 

26.5803 

1.0147 

8.6890 

1.8162 

0.7548 

0.8261 

0.7654 

18 

23 

26.9008 

0.9749 

8.1186 

1.7844 

0.7401 

0.8343 

0.7760 

18 

24 

27.2229 

1.0504 

7.5442 

1.7538 

0.7260 

0.8421 

0.7862 

18 

25 

27.5470 

1.2453 

6.9658 

1.7241 

0.7125 

0.8496 

0.7960 

18 

26 

23.2048 

2.1015 

5.7998 

1.6678 

0.6870 

0.8635 

0.8143 

18 

27 

28.6372 

3.3883 

5.0450 

1.6333 

0.6715 

0.8719 

0.8254 

18 

28 

28.7292 

4.0255 

4.8844 

1.6261 

0.6683 

0.8736 

0.8276 

19 

1 

30 .2546 

4.6668 

5.7945 

1.7407 

0.7201 

0.8454 

0.7905 

19 

2 

30.2529 

5.2955 

5.7954 

1.7407 

0.7201 

0.8454 

0.7905 

19 

3 

30.2546 

5. 9242 

5.7945 

1.7407 

0.7201 

0.8454 

0.7905 

19 

4 

30.2261 

6.5533 

5.8428 

1.7425 

0.7209 

0.8450 

0.7899 

19 

5 

30.1303 

7.1844 

5.9950 

1.7480 

0.7234 

0.8436 

0.7881 
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8.4 Output table III . - Output table IE gives the values of various variables along 
streamlines (constant ICP) selected by the input. X-SL(DIM), Y-SL(DIM), and 
Z-SL(DIM) are values of the X,Y,Z coordinates, respectively, along the streamline. The 
next four headings are the same as defined for output table II (section 8.3). The last two 
headings are the lengths (same dimension as SP(I)) of the streamline ICP computed in 
two ways as follows: 


1=1 

S-I(DIM) = X I (AX) 2 + (AY) 2 + (AZ) 2 1 

1=2 L J 


1/2 


where AX = X(I) - X(I - 1), etc. 
and 


S-II(DIM) = 



d(p 

q 


(8.4.1) 


(8.4.2) 


A sample page of output table in resulting from the sample input in section 7.4 follows: 


OUTPUT TABLE NO. Ill (COORDINATE POINTS ALONG SELECTED STREAMLINES) 


ICP 

I 

X-SL (DIM) 

Y-SL(DIM) 

Z-SL(DIM) 

Q/Q-UP 

MACH NO 

RO/RO-UP 

P/P-UP 

S-ICDIM) 

S-IICDIM) 

10 

IS 

26.5803 

9.5763 

8.6390 

1.8162 

0.7548 

0.8261 

0.7654 

27.1785 

27.2304 

10 

19 

27.9671 

9.5391 

9.4967 

1.8845 

0.7866 

0.8082 

0.7422 

28.7838 

28.8444 

10 

20 

29.2937 

9.5172 

10.3293 

1.9326 

0.8093 

0.7954 

0.7257 

30.3502 

30.4165 

10 

21 

30.5765 

9.5046 

11.1714 

1.9683 

0.8265 

0.7855 

0.7132 

31.8848 

31.9545 

10 

22 

31.8326 

9.4984 

12.0154 

1.9918 

0.8375 

0.7792 

0.7052 

33.3981 

33.4695 

10 

23 

33.0760 

9.4972 

12.8591 

2.0000 

0.8414 

0.7770 

0.7024 

34.9008 

34.9725 

10 

24 

34.3162 

9.4984 

13.7032 

2.0000 

0.8414 

0.7770 

0.7024 

36.4009 

36.4725 

10 

25 

35.5565 

9.4995 

14.5477 

2.0000 

0.8414 

0.7770 

0.7024 

37.9015 

37-9725 

10 

26 

36.7971 

9.5000 

15.3922 

2.0000 

0.8414 

0.7770 

0.7024 

39.4022 

39.^725 

10 

27 

38.0374 

9.5000 

16.2365 

2.0000 

0.8414 

0.7770 

0.7024 

40.9025 

40.9725 

10 

28 

39.2769 

9.5001 

17.0808 

2.0000 

0.8414 

0.7770 

0 .7024 

42.4023 

42.4725 

10 

29 

40.5163 

9.5003 

17.9253 

2.0000 

0.8414 

0.7770 

0.7024 

43.9023 

43.9725 

11 

1 

0.0000 

10.3C04 

6.0910 

1.0000 

0.4000 

1 . 0000 

1.0000 

0.0000 

0 . 00 0 0 

11 

2 

1.5000 

10.3002 

6.0909 

1.0000 

0.4000 

1 . 0000 

1.000 0 

1.5000 

1.5000 

11 

3 

3.0010 

10.2995 

6 .0906 

1.0000 

0.4000 

1.0000 

1.0000 

3.0010 

3.0000 

11 

4 

4.5015 

10.2935 

6.0902 

1 .0000 

0.4000 

1.0000 

1.0 0 0 0 

4.5015 

4.5000 

11 

5 

6.0020 

10.2964 

6 . 0894 

1 .0000 

0.4000 

1.0000 

1.000 0 

6.0020 

6.0000 

11 

6 

7.5024 

10.2921 

6 . 0878 

1 .0000 

0.4000 

1.0000 

1 .00 00 

7.5024 

7.5000 

11 

7 

9.0025 

10.2823 

6.0851 

1.0000 

0.4000 

1.000 0 

1 . 0 0 00 

9.0026 

9.0000 

11 

8 

10.5069 

10.2566 

6.0810 

1 .0132 

0.4055 

0.9979 

0 . 9970 

10.5072 

10.5049 

11 

9 

12.0288 

10.2086 

6 . 0801 

1 .0490 

0.4203 

0.9920 

0.9888 

12.0298 

12.0233 

1 X 

10 

13.5837 

10.1430 

6.0961 

1.1007 

0.4413 

0.9832 

0.9765 

13.5862 

13.5855 

11 

11 

15.1835 

10.0662 

6.1518 

1.1619 

0.4674 

0 . 9722 

0.9613 

15.1883 

15.1831 

11 

12 

16.8276 

9.9757 

6.2757 

1 .2365 

0.4988 

0.9582 

0.9420 

16.8401 

16.3330 

11 

13 

18.4936 

9.8688 

6.4955 

1 .3296 

0.5385 

0 . 9397 

0 . 9166 

18.5239 

18.5215 

11 

14 

20.1495 

9.7535 

6.8298 

1 .4361 

0.5845 

0 . 9172 

0.8860 

20.2171 

20.2182 

il 

15 

21.7684 

9.6423 

7.2840 

1.5493 

0.6341 

0.8917 

0.8518 

21.9022 

21.9112 

11 

76 

23.3319 

9.5457 

7.8511 

1 .6613 

0.6841 

0.8651 

0.8164 

23.5681 

23.5879 

11 

±7 

24.8302 

9.4692 

8.5152 

1 .7640 

0.7307 

0.8395 

0.7828 

25.2088 

25.2402 

11 

18 

26.2611 

9.4142 

9.2552 

1 .8491 

0.7700 

0.8176 

0.7543 

26.8207 

26 .8627 

11 

19 

27.6287 

9.3790 

10.0476 

1 .9081 

0.7977 

0.8019 

0.7342 

28.4017 

28.4522 

11 

20 

28.9431 

9.3581 

10.8699 

1 .9467 

0 .8159 

0.7915 

0.7209 

29.9523 

30 .0088 

11 

21 

30.2195 

9.3459 

11.7059 

1.9754 

0.8296 

0.7837 

0.7109 

31.4731 

31.5387 

11 

22 

31.4729 

9.3399 

12.5468 

1.9936 

0.8383 

0.7788 

0.7046 

32.9875 

33.0504 

11 

23 

32.7158 

9.3386 

13.3892 

2.0000 

0.8414 

0.7770 

0.7024 

34.4889 

34.5529 

11 

24 

33.9559 

9.3397 

14.2327 

2.0000 

0.8414 

0.7770 

0.7024 

35 . 9887 

36 . 0529 

11 

25 

35.1962 

9.3408 

15.0770 

2.0000 

0.8414 

0.7770 

0.7024 

37.4891 

37.5529 

11 

26 

36.4367 

9.3412 

15.9216 

2.0000 

0.8414 

0.7770 

0.7024 

38.9398 

39.0529 

11 

27 

37.6768 

9.3412 

16.7659 

2.0000 

0.8414 

0.7770 

0.7024 

40.4901 

40.5529 

11 

28 

38.9162 

9.3412 

17.6102 

2.0000 

0.8414 

0.7770 

0.7024 

41 . 9393 

42.0529 

11 

29 

40.1557 

9 . 3414 

18.4550 

2.0000 

0.8414 

0.7770 

0.7024 

43.4898 

43.5529 

12 

1 

0 . 0 0 0 0 

9.2730 

7.6320 

1 .0000 

0.4900 

1.0000 

1.0000 

0.0000 

0.0000 

12 

2 

1.5000 

9.2727 

7.6818 

1.0000 

0.4000 

1.0000 

1.0000 

1.5000 

1.5000 

12 

3 

3.0009 

9.2720 

7.6814 

1.0000 

o.^.ooo 

1.0000 

1.0000 

3.0009 

3.0000 

12 

4 

4.5014 

9.2707 

7.6803 

1.0000 

0.4000 

1.0000 

1.0000 

4.5014 

4.5000 

12 

5 

6,^019 

9.2682 

7.6793 

1.0000 

0.4000 

1 .0000 

1.0000 

6 .0019 

6.0000 
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8.5 Output to tape or disk for three-dimensional graphics . - Output to tape or disk 
for three-dimensional graphics occurs at the end of subroutine PUTOUT and is described 
in section 6.3. A three-dimensional plot resulting from the sample input in section 7.4 
follows. Running time on an IBM 370/3033 was 6.50 min. 



9.0 NUMERICAL EXAMPLES 

Five numerical examples of ducts are presented. For each example, the upstream 
boundary configuration and associated grid are given together with the prescribed velocity 
distribution on the lateral surface and a number of key input parameters. The results are 
presented by three-dimensional graphs. The first example is a completely general three- 
dimensional nozzle with a nonsymmetrical upstream boundary configuration and rapid 
acceleration of the flow with no deceleration along the surface streamlines. The second 
example is an accelerating elbow with the same upstream boundary configuration and 
again no deceleration along the surface streamlines. The third example is an accelerating 
S-duct with an elliptical upstream boundary configuration. The fourth example is a 
rapidly decelerating elbow with a circular upstream boundary and an unusually sharp 
turning angle. This solution, like the others, can be reversed to give, in this case, a 
rapidly accelerating elbow with no deceleration along the surface streamlines. Of special 
interest in this example is the pronounced initial turning of the inner wall in a direction 
opposite to that of the elbow itself. This phenomenon has also been observed (ref. 3) in 
designs of two-dimensional ducts. The last example is a preliminary design of a side-inlet 
duct such as might be used with various types of turbomachinery. The solution has planar 
symmetry (with a small amount of overlap in one region), and for the reverse flow case, is 
an accelerating flow into a circular annulus with no deceleration anywhere along the duct 
walls. 
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9.1 Numerical example I . - Straight, three-dimensional nozzle with rapid acceleration (no 
deceleration along streamlines) 





I L - 31, I 2 - 37, I 3 ■ 31, I 4 - 37, NI - 67, q D = 2.0, AS - 0.5 


Option ISYM = 1 Accuracy of finite-difference solution (EPSR) * 0.000005 

Option IVEL = 1 Overrelaxation factor (ORELAX) =1.35 

Major iterations (ITER) = 4 Exit-area error (ERRAR) = 0.0005 

Coefficient to average x, y, z(CAVP) = 0.5 Running time (370/3033), min = 55.69 

Upstream Mach number (AMU) = 0,4 DEL-P-01 = 0.3 

Ratio of specific heats (GAM) = 1,4 DEL-P-12 = 0.2 

J location of primary streamline (JX) = 9 DEL-P-23 - 0.3 

K location of primary streamline (KX) = 12 Lxation of principal streamline (XP) » 0.9 

Number of subdivisions between 
adjacent input values of I (NSD) = 1 


56 



numerical EXAMPLE 


r*~ 

CO 



9.2 Numerical example II . - General case of three-dimensional accelerating elbow (no 
deceleration along streamlines) 
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Accuracy of finite-difference solution (EPSR) 

= 0.000005 

Overrelaxation factor (ORELAX) 

■ 1.30 

Exit-area error (ERRAR) 

- - 0.0019 

Running time (370/3033), min 

- 103.14 

DEL-P-01 

* 0.3 

DEL-P-12 

■ 0.2 

DEL-P-23 

- 0.3 

Location of principal streamline (XP) 

• 0.9 







9.3 Numerical example III . - Accelerating S-duct with elliptical upstream boundary (no 
deceleration along streamlines) 
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NUMERICAL EXAMPLE III - UPSTREAM BOUNDARY 
AND ASSOCIATED GRID 
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Option ISYM 
Option IVEL 
Major iterations (ITER) 

Coefficient to average x, y, z (CAVP) 
Upstream Mach number (AMU) 

Ratio of specific heats (GAM) 

J location of primary streamline (JX) 
K location of primary streamline (KX) 
Number of subdivisions between 
adjacent input values of I (NSD) 


1 

3 

10 

0.2 

0.3 

1.4 

4 
8 


Accuracy of finite-difference solution (EPSR) = 0. 000005 
Overrelaxation factor (ORELAX) = l. 30 

Exit-area error (ERRAR) = 0.0005 

Running time (370/3033), min = 29.64 

DEL-P-01 = * 

DEL-P-12 = * 

DEL-P-23 = * 

Location of principal streamline (XP) • * 


1 


* Not applicable. 
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9.4 Numerical example IV . - Decelerating elbow with sharp turn and circular upstream 
boundary (no deceleration for reversed flow) 
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INPUT FOR NUMERICAL EXAMPLE IV 



Option ISYM = 1 

Option IVEL ■ 2 

Major iterations (ITER) = 4 

Coefficient to average x, y, z(CAVP) = 0.0 
Upstream Mach number (AMU) * 0,6 

Ratio of specific heats (GAM) * 1,4 

J location of primary streamline (JX) = 10 
K location of primary streamline <KX) * 10 
Number of subdivisions between 
adjacent input values of I (NSD) = I 


Accuracy of finite-difference solution(EPSR) = 0.000005 
Overrelaxation factor (ORELAX) = 1.30 

Exit-area error (ERRAR) = 0.0024 

Running time (370/3033), min = 121.88 

DEL-P-01 = * 

DEL-P-12 - * 

DEL-P-23 - * 

Location of principal streamline (XP) - * 


Not applicable. 





9.5 Numerical example V . - Planar symmetry solution for side inlet (in reversed-flow 
case; no deceleration for reversed-flow direction) 





INPUT FOR NUMERICAL EXAMPLE V 



1 I I + 1 I x I 2 I3 I4 MI 


Ij ■ 21, I 2 ■ 31, I 3 - 11, I 4 = 31, NI = 60, q D - 0. 25, AS = 2.0 

Option ISYM = 2 Accuracy of finite-difference solution(EPSR) = 0.000005 

Option IVEL ■ 2 Overrelaxation factor (ORELAX) = 1,30 

Major iterations (ITER) - 6 Exit-area error (ERRAR) - 0.0110 

Coefficient to average x, y, z(CAVP) ■ o.O Running time (370/3033), min = 33.05 

Upstream Mach number (AMU) = 0.6 DEL-P-01 = * 

Ratio of specific heats (GAM) = 1.4 DEL-P-12 = * 

J location of primary streamline (JX) = 9 DEL-P-23 = * 

K Ixation of primary streamline (KX) ■ 15 Location of principal streamline (XP) = * 

Number of subdivisions between 
adjacent input values of I (NSD) * 1 

Not applicable. 
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numerical example V 
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10.0 CONCLUDING REMARKS 

The general design method for three-dimensional, potential flow developed in part I 
of this report (ref. 1) is herein applied to the design of simple, unbranched ducts. A com- 
puter program, DIN3D1, is developed and five numerical examples are presented, 
including a nozzle, two elbows, an S-duct, and the preliminary design of a side inlet for 
turbomachines. The two major inputs to the program are the upstream boundary config- 
uration and the lateral velocity distribution on the duct wall. As a result of these inputs, 
boundary conditions of the problem are overprescribed and the problem is ill posed. How- 
ever, it appears that there are degrees of "compatibility" between the two major inputs 
and that for reasonably compatible inputs satisfactory, reliable solutions can be obtained. 
By not prescribing the shape of the upstream boundary, the problem presumably becomes 
well posed, but it is not clear how to carry out a practical design method under this cir- 
cumstance. Nor does it appear desirable, because the designer usually needs to retain 
control over the upstream (or downstream) boundary configuration. 

The problem is further complicated by the fact that, unlike the two-dimensional 
case, and irrespective of the upstream boundary shape, some prescribed lateral velocity 
distributions do not have proper solutions (appendix C). 

The input data for an example solution together with example output tables and a 
three-dimensional plot of the solution are given in sections 7.4 and 8.1 to 8.5, respectively. 
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APPENDIX A 


"EQUILIBRIUM" VELOCITY DISTRIBUTIONS FOR 
INPUT OPTION IVEL EQUAL TO 2 OR 3 

The "equilibrium" velocity distributions for input option IVEL equal to 2 or 3 refer 
to the velocity distributions around the contours of the potential surfaces; the velocity 
distribution along the principal streamline (section 7.2) is not affected. Variation in the 
velocity distribution around the contour (e.g., DQ in second fig. of section 7.2) causes the 
duct to bend and may be looked upon as the duct "loading." 

Consider potential flow in an infinitely long duct with constant loading. Such a duct 
will turn an infinite number of degrees, and the duct cross section will be constant. Under 
these circumstances, the potential surfaces are flat planes, and the "equilibrium" velocity 
distribution normal to the planes is a free vortex 

qr = constant (Al) 

where the radius r is measured from the axis about which the duct bends. 

Such an equilibrium duct shape halfway between +°° can be considered to lie on 
the Y,Z plane corresponding to the upstream boundary. Thus, 



For IVEL = 2, the axis of the bend is a line of constant Z; and for IVEL = 3, the axis is a 
line of constant Y, as shown. 

For IVEL equal to 2 or 3, the "equilibrium" shape is assumed to be the input shape 
of the upstream boundary and the lateral velocity distribution corresponds to the "equi- 
librium" velocity based on that shape. (Other shapes could be used, but these would entail 
additional input and probably would not achieve the same degree of compatibility 
(section 4.1) between the prescribed upstream boundary shape and the prescribed lateral 
velocity distribution.) 



(A2) 


For IVEL = 2, equation (Al) gives 

’ r "Vp -( < *P + A «.mp) r o 

where (preceding figure) the subscript P refers to the principal streamline (at which qp 
is the input value of QP(I)) and subscript o refers to the maximum (outer) radius where 
q is equal to qp + Aq amp and Aq amp is the input value ( negative) of DQAMP(I) at 
potential surface I. (Note that DQAMP(I) for IVEL options 2 and 3, as opposed to 
option 1, is the difference between the minimum velocity, which occurs at r 0 in the 
figure, and the input velocity QP(I) of the principal streamline.) 


In the figure, the radius r is related to the Z coordinate by 


r = r„ + r 
P l o 



(A3) 


and 


r 

o 



- Z 

o 


From equations (A2) and (A4), 

Vp-(% +A ’amp)[ r P + ( Z p- Z o)] 

or 


r P = 


q p + Aq 

r amp 

- Aq 

amp 



and from equations (A2), (A3), and (AS), 


(A4) 


(A5) 
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(A6) 


„ Vp l? 

•' r ‘ , f 2 ?- 2 \ 

q P + A<5 amp \ Z ?- Z o) 

Thus, at each potential surface I, for IVEL = 2, equation (A6) gives the lateral distri- 
bution of velocity q as a function of the coordinate Z around the upstream boundary 
contour. 

As shown in the preceding figure, rp (at which radius q = qp) can be greater 
than r -mir, (at which radius q = q max )> but r„ must be reasonably greater than 
r 0 (at which radius q = q-min )- In subroutine VELD, for IVEL = 2, if 



r - r . "Z . - Z 
o mm nun o 


< 0.5 


the solution is stopped. 

In a similar fashion, for IVEL = 3, 


q P r P 


Aq a 


mp_ 


Y - Y 
x p 


q P + Aq amp V Y P- 


and the solution is stopped if 


V r P y p- y o 


r - r . Y . - Y 

o nun nun o 


< 0.5 


(A7) 


(A8) 


(A9) 
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APPENDIX B 


CONDITION FOR NORMALITY OF UNIT VECTORS e 2 and e 3 
WITH UNIT VECTOR e x 

Consider the case in which the unit vectors e]_ and e 2 are not normal and find the 
direction cosines for a third unit vector e 2x , which lies in the plane of ei and e 2 and is 
normal to ei. Thus, 


e 2x 


Because the three vectors are coplanar, they are related to 

v * V i * ¥2 

from which 

cos cos cos 

cos B^ 

C0S Y 2x = \ C0S Y 1 + k 2 C0S Y 2 


cos B„ 


cos B„ 


(Bl) 


(B2) 
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where cos a 2 X ,...» cos Y2 are the direction cosines of e2x? e l> and e 2 and and k 2 
are constants. 

The constants kj and k 2 are determined from equation (B I) as follows: 

? 2x • '! -°- k I * k 2 C0S9 

and 

F„ . F„ = cos (9 - 90°) = sin 9 = k, cos 9 + 1c 
2x 2 12 

from which 


^ 1 tan 9 

and 


(B3) 


^2 sin 9 

Thus, the direction cosines for F 2 X , which is F 2 adjusted to satisfy the normality 
condition, are known from equations (B2) to (B4). 


(B4) 


In a similar fashion, the adjusted direction cosines for eg are given by 


cos o^ x = k^ cos 
cos 

cos Y 3x « \ cos y 1 


cos . k. 



COS 


*3 


cos 


cos Y 0 


(B5) 


where k^ and k 2 are given by equations (B3) and (B4), respectively. 
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APPENDIX C 


PROBLEMS AND LIMITATIONS OF THE THREE-DIMENSIONAL DESIGN METHOD 

The ill-posed nature of the three-dimensional duct design method when both the 
upstream boundary configuration and the lateral velocity distribution are prescribed is 
discussed in section 4.0. This ill-posed nature negates a proper solution. However, for 
relatively "compatible" upstream boundary configurations and lateral velocity distribu- 
tions, reasonable solutions are forced by introduction of the six constraints in section 4.2 
and by limiting the number of major iterations (ITER) to a range over which the solution is 
converging, as evidenced by decreasing maximum residuals 3t. This appendix considers 
other problems and limitations of the design method, and for this purpose it is assumed 
that a method exists for assuring an absolutely compatible upstream boundary configu- 
ration (if such exists) for a given prescribed, lateral- velocity distribution. Also, without 
destroying the generality of the discussion, it is convenient to assume planar symmetry. 

For the classical, two-dimensional , duct design problem (ref. 3) a solution exists for 
every prescribed, piecewise- continuous, velocity distribution along the duct walls. For 
the three-dimensional problem this universal existence does not appear to be the case. 

For example, consider a straight duct (which, if the velocity is not constant, implies two 
planes of symmetry at right angle, with the duct centerline along the intersection). 
Presumably, if the duct is straight, the upstream boundary configuration has biplanar 
symmetry, but is otherwise general. Thus, 



where P is the decimal fraction of distance around the contour. In the upstream and 
downstream regions (section 7.2), the lateral velocity q, expressed as a ratio of the 
upstream velocity, is 1.0 along every boundary streamline. Elsewhere, let the prescribed 
velocity be 1.0 along the streamlines through a and c (P = 0.0 and 0.5, respectively), and 
let the velocity decrease along the contour in an arbitrary fashion, but with biplanar sym- 
metry, to a finite value approaching zero for the streamline through b (P = 0.25). If this 
distribution of velocity with P is maintained along the boundary streamlines over a large 
range of the velocity potential Acp, then from equation (3.0.1) 


80 



Acp = q (As) = q,(As) 
3-jC a,C D D 


(Cl) 


so that over the range Acp 



(C2) 


from which the length (As)^ of streamline b becomes many times larger than the 
streamline length (As) a>c , and no solution (i.e., shape of flow field) appears likely. It 
might be argued that the large (As)fo could be accommodated by a rapid outward fanning 
of streamline b, but the pressure gradients associated with the velocity distribution 
preclude this. (The rapid outward fanning of streamline b would approach a two- 
dimensional configuration in which streamlines a and c come together, lie on the duct 
centerline, and have a velocity distribution that adjusts to the prescribed velocity of 
boundary streamline b. In the three-dimensional case, however, the velocity distributions 
along streamlines a and c are prescribed and thus cannot adjust. This inability to adjust 
is probably the center of the problem.) 


Finally, for the example just discussed the velocity q along the straight streamline 
on the centerline of the duct can be no higher than 1.0 (which is the highest velocity on 
the boundary streamlines in this example) and will be less than 1.0 where influenced by 
velocities less than 1.0 on the boundary. The lengths of the streamlines between the 
upstream and downstream potential surfaces are given by 



Thus, As for the centerline streamline with velocities less than 1.0 is longer than As for 
the boundary streamlines a and c, which have a constant prescribed velocity of 1.0. 
However, in contradiction, the centerline streamline must be shorter than the boundary 
streamlines a and c, because it is straight and normal to the upstream and downstream 
potential surfaces, which are flat and parallel. It is concluded that for three-dimensional 
design problems not every prescribed velocity distribution has a proper solution. 

For velocity distributions without proper solutions program DIN3D1, using the 
constraints in section 4.2 and limiting the number of major iterations (ITER), forces a 
"reasonable" solution. A measure of this reasonableness is the difference in streamline 
lengths S-I and S-II (output table III), which lengths should be equal. 

Another problem area in the application of program DIN3D1 occurs when certain 
characteristics in the shape of the downstream boundary configuration are desired. 

(There is, of course, no way to achieve a precise shape, because the downstream config- 
uration is dictated by the upstream configuration and the prescribed lateral velocity dis- 
tribution.) For example, consider a straight duct with a transition section in which the 
duct cross section changes from a circular upstream shape to an elliptical downstream 
shape of the same area. A normal design procedure, based on one -dimensional consider- 
ations, would keep the duct area constant and employ a linear variation in the fineness 
ratio of the elliptical cross section starting from 1.0 for the circle and ending with the 
desired value for the downstream shape. Here, to avoid large losses, the designer's 
objective is to keep the velocity on the duct wall constant (i.e., q = 1.0 along the boundary 
streamlines); and provided that the transition length is not too short, this objective should 
be nearly achieved. Thus for this type of three-dimensional design problem very large 
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changes in the duct cross section occur for very small changes in the prescribed lateral 
velocity distribution. It is not easy to determine just how large and where these small 
changes in velocity should be. Furthermore because this type of problem is so sensitive to 
small changes in lateral velocity distribution, the downstream boundary shape is also 
sensitive to the necessarily approximate methods used in the forced, finite-difference 
solution of the governing differential equation (3.3.1). That is, the downstream shape 
varies appreciably with number of iterations (ITER), with upstream grid size and arrange- 
ment (section 7.1), with various damping coefficients (sections 4.6 and 5.2.2), and perhaps 
with such lesser iterations as ICONX (section 5.1.1) and NTRY (section 7.3). 

In summary, it is not easy to control the downstream boundary shape of the duct by 
the prescribed lateral velocity distribution; although it is relatively easy to control the 
streamwise shape of the duct by this means. However, substantial differences in the 
downstream shape need not imply significantly different lateral velocity distributions, 
provided that the downstream areas are equal in size. Finally, in those cases where the 
downstream shape is important, the solution becomes more sensitive to the lateral veloc- 
ity distribution, if the damping coefficient CAVN is reduced to 0.0, or at most is not 
greater than 0.1. 
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